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Abstract 

Background: The past thirty-five years have seen an intense search for the molecular 
mechanisms underlying calcium-induced calcium-release (CICR) in cardiac myocytes, 
with voltage clamp (VC) studies being the leading tool employed. Several VC 
protocols including lowering of extracellular calcium to affect Co 2+ loading of the 
sarcoplasmic reticulum (SR), and administration of blockers caffeine and thapsigargin 
have been utilized to probe the phenomena surrounding SR Co 2+ release. Here, we 
develop a deterministic mathematical model of a rat ventricular myocyte under VC 
conditions, to better understand mechanisms underlying the response of an isolated 
cell to calcium perturbation. Motivation for the study was to pinpoint key control 
variables influencing CICR and examine the role of CICR in the context of a 
physiological control system regulating cytosolic Co 2+ concentration ([Ca 2+ ] myo ). 

Methods: The cell model consists of an electrical-equivalent model for the cell 
membrane and a fluid-compartment model describing the flux of ionic species 
between the extracellular and several intracellular compartments (cell cytosol, SR and 
the dyadic coupling unit (DCU), in which resides the mechanistic basis of CICR). The 
DCU is described as a controller-actuator mechanism, internally stabilized by negative 
feedback control of the unit's two diametrically-opposed Co 2+ channels (trigger- 
channel and release-channel). It releases Co 2+ flux into the cyto-plasm and is in 
turn enclosed within a negative feedback loop involving the SERCA pump, 
regulating [Co 2+ ] myo . 

Results: Our model reproduces measured VC data published by several laboratories, 
and generates graded Co 2+ 
environment where [Ca 

the DCU elements in this process, particularly the role of the ryanodine receptor in 
controlling SR Co 2+ release, its activation by trigger Co 2+ , and its refractory 
characteristics mediated by the luminal SR Co 2+ sensor. Proper functioning of the 
DCU, sodium-calcium exchangers and SERCA pump are important in achieving 
negative feedback control and hence Co 2+ homeostasis. 

Conclusions: We examine the role of the above Co 2+ regulating mechanisms in 
handling various types of induced disturbances in Co 2+ levels by quantifying cellular 
Co 2+ balance. Our model provides biophysically-based explanations of phenomena 
associated with CICR generating useful and testable hypotheses. 



release at high Co gain in a homeostatically-controlled 
myo is precisely regulated. We elucidate the importance of 
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Background 

Contraction of cardiac muscle is triggered by a transient rise in intracellular Ca 2+ con- 
centration [Ca 2+ ] myo . Sarcolemmal (SL) membrane depolarization triggers Ca 2+ influx 
from the extracellular medium by opening dihydropyridine (DHP)-sensitive L-type 
Ca 2+ channels. Following diffusion across a small sub-membrane dyadic space, this 
influx activates ryanodine receptors (RyRs) controlling ryanodine-sensitive Ca 2+ release 
channels in the junctional portion of the sarcoplasmic reticulum (jSR). Fabiato and 
Fabiato [1] named the process calcium-induced calcium release (OCR). Ca 2+ subse- 
quently diffuses from the dyadic space into the cytosol Ultimately, intracellular Ca 2+ 
concentration [Ca 2+ ] myo is returned to resting levels by combination of: (a) Ca 2+ buf- 
fering in the dyadic space and cytosol; (b) sequestration of Ca 2+ by sarcoplasmic/endo- 
plasmic reticulum Qz 2+ -ATPase (SERCA)-type calcium pumps lining the longitudinal 
portion of the sarcoplasmic reticulum (LSR); and (c) Ca 2+ extrusion from the cytosol 
by Na + /Ca 2+ exchangers and Qz 2+ -ATPase pumps on the sarcolemmal membrane. 

CICR in cardiac muscle exhibits both graded behavior and a high gain. Graded beha- 
vior refers to the obser-vation that SR Ca 2+ release is proportional to the influx of trigger 
Ca 2+ [2], whereas high gain indicates that the SL trigger current elicits a high SR Ca 2+ 
release flux. Graded Ca 2+ release with high gain is somewhat paradoxical according to 
Stern [3], in that the positive feedback inherent in such high-gain systems tend to pro- 
duce regenerative, nearly all-or-none release rather than graded release. Several determi- 
nistic models have been developed to explain excitation-contraction (E-C) coupling 
[4,5], but none of them can explain the mechanism of graded release at high gain over a 
wide range of values for sarcolemmal Ca 2+ current. Stern [3] proposed that such a gra- 
dation paradox might be explained if the stimulus for Ca 2+ release by RyRs were actually 
the local nanodomains of [Ca 2+ ] generated by nearby L-type channels, rather than the 
global cytosolic [Ca 2+ ] myo . According to this hypothesis, graded control of macroscopic 
SR Ca 2+ release can be achieved by graded statistical recruitment of individual, autono- 
mous, all-or-none stochastic release events [6]. In these studies, a distributed differential 
model of high order that included dynamic interactions between large numbers of indi- 
vidual channels was used to demonstrate this concept. However, rather large amounts of 
computation time are required with distributed stochastic models of this type. Addi- 
tional models have sought to characterize the Ca 2+ release complex, including several 
[7-9] based on the stochastic release process adopted by Stern et al. These statistical 
models have solved the graded release problem, however, they too are complicated and 
computationally very expensive. Other models based on the simplified local control 
model of CICR developed by Hinch et al. [10] sought to adopt a lower order description 
of the E-C coupling process [11,12] by making an approximation of rapid equilibrium in 
the dyadic space. The latency from onset of Ca 2+ entry via the I Ca ,L channel to triggered 
SR Ca 2+ release is known to increase with decrease in the magnitude oil Ca ,L,TT [13], the 
modeling of which is made possible by considering Ca 2+ diffusion in the dyadic medium. 
These models [11,12] also approximate the SR as a single volume compartment with no 
distinction between junctional versus the longitudinal (network) SR compartments. 
However, recent work [14] points towards the important role of the Ca 2+ refilling rate 
from the network to junctional SR in controlling RyR release termination via the luminal 
sensor. Shiferaw et al. [15] developed a computationally tractable model of Ca 2+ cycling 
to represent the release of calcium from the SR as a sum of spatially localized events 
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that correspond to Ca + sparks, assuming the recruitment rate of Ca + sparks is directly 
proportional to the whole-cell I Ca ,L current. This assumption overlooks the complex cal- 
modulin mediated interaction (calcium dependent facilitation (CDF) and calcium depen- 
dent inactivation (CDI)) of the I Ca ,L channel with calcium in its vicinity. It also demands 
a large amount of computation. 

Numerically, distributed as well as statistical models tend to be computationally 
expensive due to the in-herent repetition involved in the computation. In a spatially 
distributed model, simultaneous solution for dynamics in identical compartments dis- 
tributed in space would amount to a large computational cost. In statistical models 
inference is drawn based on multiple runs of identical events which translate into a 
prolonged simulation time. Hence, these models are cumbersome to implement, parti- 
cularly in larger multiple-cell simulations. Consequently, we consider a deterministic 
approach to the characterization of CICR. Specifically, we develop a lumped model of 
the Ca 2+ release complex that includes: (a) a sub-sarcolemmal dyadic cleft space separ- 
ating the SL and jSR membranes; (b) a single DHP-sensitive Ca 2+ channel on the SL 
membrane; and (c) a single equivalent Ry-sensitive channel arranged symmetrically on 
the opposing jSR membrane that represents the output of a local cluster of Ry-sensitive 
channels facing the DHP-sensitive channel. Based on morphological data compiled by 
Bers [16], we further assume that each ventricular cell contains 10,000 of these dyadic 
Ca 2+ release units, and that they are associated with the fraction of the SL membrane 
that is coupled with the jSR. That is, we partition the sarcolemma into free and dyadi- 
cally coupled SL membrane, and associate each with a different fluid compartment: the 
cell cytosolic medium in the case of the free SL membrane, and the dyadic cleft space 
medium in the case of the dyadic-coupled fraction. In a sense, we build on Stern's [3] 
local domain concept by considering the aforementioned local nanodomains identical, 
but focusing on the nonlinear dynamics of the two different types of Ca 2+ channel in 
the dyadic coupling unit. Our deterministic model although is very descriptive, is com- 
putationally tractable and has a run time of 21 sec (including recording of 73 variables 
of type double on a data file) for 1 cycle of 4 Hz voltage clamp stimulation. 

Methods 

Experimental Methods 

Rat ventricular myocytes were prepared from 200-300 g male Sprague Dawley rats by 
dissociation with collagenase, as previously described [17]. All experiments were per- 
formed under conventional whole cell recording conditions with a List EPC-7 patch 
clamp, recording fluorescence from nearly the entire cell, as described by Fan and 
Palade [17]. Recordings from an individual cell were rarely extended beyond 10 min in 
order to reduce as much as possible both escape of dye from the cell and Ca 2+ current 
rundown. External solution in the bath was normal Tyrode (1 mM Ca 2+ ) with Cs + sub- 
stituted for IC for purposes of blocking inward rectifier IC currents. The internal solu- 
tion in the pipette contained Cs aspartate supplemented with 20 mM CsCl, 3 mM 
Na 2 ATP , 3.5 mM MgCl 2 and 5 mM HEPES. Holding potential used was -40 mV. 

Computational Aspects 

All simulations and analysis were performed on a 2.8 GHz Intel® Core™2 Duo CPU- 
based computer using Microsoft Windows XP operating system. To find the 
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parameters involved in the 6 state Markovian model for Ic a ,u a non-linear least- 
squares method [18] was used to obtain the solution of the system of non-linear ordin- 
ary differential equations. Specifically, we have employed an algorithm given by Lau 
[19]. The numerical integration scheme used to solve the full set of forty two lst-order 
differential equations describing the dynamic model was the Merson-modified Runge- 
Kutta 4th-order method [20,21] with a conservative fixed time step, chosen small 
enough to allow the local truncation error to be of fourth order. The explicit finite dif- 
ference scheme was used to numerically solve the Laplacian equations of Ca 2+ diffu- 
sion in the cleft space. Detailed numerical methods are similar to those presented by 
Smith et al. [22]. The results were visualized using Matlab by Mathworks and Origin 
by Microcal Software. 

Model Development 

Our objective was to develop a model of the rat ventricular cell which could be used to 
explain Ca 2+ signaling at the nanoscale level of the dyad and integrate the contribu- 
tions of many dyads to produce a Ca 2+ transient and continuous Ca 2+ balance at the 
whole-cell level. Therefore, we start with a broad discussion of the elements of the 
DCU and its Ca 2+ supply (the jSR), and continue with a progressively more detailed 
description of the whole cell model. It is important to note that all Ca 2+ concentra- 
tions discussed in the model pertain to unbound Ca 2+ unless specified. 

Membrane Classification 

We assume that a continuous membrane barrier exists between the cytoplasm and the 
external bathing medium (Figure 1A; Figure 2), which consists of two components: a 
surface sarcolemma (SL) {M FreeS L in Figure 3) free of any sub-membrane contact with 
the junctional sarcoplasmic reticulum (jSR) and the remnant membrane {M Junctiona i S L 
in Figure 3) that does make contact with the jSR via a dyadic space (nanodomain) (Fig- 
ure 1A). These membrane components have the same basic plasma membrane, but dif- 
fer in content with regard to total membrane surface area, type and distribution of 
transmembrane ion channels, ATPase pumps and exchangers, as well as their func- 
tional coupling with a dyadic space. Ultrastructural information from several cardiac 
preparations including the rat ventricular cell has been compiled by Bers [16], which 
can be used to estimate the percentage of the cell membrane in contact with a dyadic 
space, for either the free surface plasmalemma or for the transverse tubule (TT) which 
brings the extracellular medium to the plasma membrane of the dyadic coupling unit. 
Thus, the bounding membrane is divided into two lumped parts (free and coupled) 
based on the existence of sub-membrane coupling to a dyadic space (Table 1). A por- 
tion of membrane could be part of a transverse tubular membrane, but if there is no 
dyadic coupling involved, that membrane would be classified as belonging to the free 
surface plasmalemma. Another portion of membrane might be part of the bounding 
outer surface of the cylindrical cell and yet have submembrane coupling to a dyadic 
space. In this case, it would be classified as belonging to the coupled category. Table 2 
gives values for volumes of the fluid compartments (shown in Figure 3) assumed for 
the rat ventricular cell, which are largely based on measured data from rat ventricular 
myocytes [16]. 
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Figure 1 Cellular fluid compartments. Figure 1: Cellular fluid compartments, (a) Model configuration 
showing dyadic space, jSR, LSR, cytoplasm and SL; (b) Inset provides a more detailed description of the 
dyadic space showing the coupling of the two types of Co 2+ channels (trigger and Co 2+ release channels) 
via the dyadic fluid medium. Only one representative dyadic coupling unit is shown, however the whole 
model contains such 10,000 identical units lumped together. 




Dyadic Space 



Channel and Exchanger Distribution 

Recent research has also shown that besides L-type Ca 2+ channels, Na + ICa + exchanger 
activity is also found predominantly in the T-tubules of rat ventricular myocytes [23]. 
Our model configuration reflects this finding in that the tubular fraction of I Ca ,L ^NaCa 
and ImuCs channels facing a unitary dyadic space are denoted as ica,L, tt (source of trig- 
ger Ca 2+ into a unitary dyad), ijsiaCa, tt and iNaCs,TT respectively (Figure 2). The free 
sarcolemmal component of these same channels are denoted as Ica,L,su L NaCa,SL an d 
In<zCs,sl respectively. 1cu,l,tt is the total current entering through L-type Qz 2+ channels 
via all the dyadic units (N dyad ). We define the total L-type current I Ca ,L as the combi- 
nation of l Ca ,L,TT and l Ca> L,SL (Le., l Ca> L = Ica,L,TT + Ica,L,sd- ^Ca,L in our model is mostly 
(90%) from the L-type Ca 2+ current in the T-tubules, since Kawai et al. [24], found L- 
type current to be highly concentrated (9-fold) in the T-tubules (Icci,l,tt) vs. the cell 
surface sarcolemma {Ic a ,L,SL.) of rat ventricular myocytes. We described the l Ca ,L chan- 
nel using a 6-state Markovian model as shown in Figure 4A. The distribution of I^aCa 
and ImuCs correspond to that of I Ca> L channel in order to ensure Ca 2+ and Na + ion 
balance. 

Since our study is focused on voltage clamp testing of Ca 2+ transients in rat ventri- 
cular cells, we assume that the majority of Na + and IC channels are blocked by either 
the holding potential used (-40 mV) or appropriate blocking agents. Thus, these chan- 
nels are not modeled, and we assume that only the dihy-dropyridine (DHP) -sensitive 
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Figure 2 Electrical equivalent circuit for the plasma membrane of a rat ventricular cell. Figure 2: C m 
77 : membrane capacitance of the junctional SL membrane coupled with the dyadic space; C mi5L : 
membrane capacitance of the uncoupled free SL membrane; Currents through the uncoupled free SL 
membrane are (a) \co,l,sl- L-type calcium current, (b) \ NOiS l- sodium current through the DHPR channel, (c) l Cs , 
si. cesium current through the DHPR channel, (d) l Na ca,sL' sodium-calcium exchanger current, (e) \ Na cs,sL' 
sodium-cesium exchanger current, (f) \pmca- calcium pump current, (g) \ NOib : background sodium current- 
Currents through the junctional SL membrane coupled with a dyadic space are (h) \ COi ljt '■ L-type calcium 
current; (i) 1^,77 : sodium current through the DHPR channel, (j) \ CsJT : cesium current through the DHPR 
channel, (k) \ Na cajr '■ sodium-calcium exchanger current, (I) \ Na csjr '■ sodium-cesium exchanger current; V G : 
potential in external medium; V,-: intracellular potential; The coupling resistance between the surface 
sarcolemma and the transverse tubules being very small is neglected in our model hence V m : common 
transmembrane potential across both uncoupled and coupled membranes. 



Ca 2+ channels, the electrogenic pumps, Na + /Ca + exchangers and Na + /Cs + pumps 
expressed in the free and/or coupled SL membranes contribute to the voltage clamp 
response. Table 3 provides values for various parameters used to model the ion trans- 
port across the sarcolemmal membrane. 

The SR Fluid Compartment 

The SR is an intracellular organelle that consists of two lumped fluid compartments 
(the jSR and LSR) that communicate (Figure 1A; Figure 3). Like the sarcolemma, the 
bounding membranes of the jSR and LSR are differentiated regarding their ionic cur- 
rent content and degree of coupling with the sarcolemma. With regard to ionic cur- 
rents, the LSR membrane has a thapsigargin-sensitive SERCA pump for pumping Ca 2+ 
into the LSR lumen against a concentration gradient. In contrast, the jSR membrane 
contains an outwardly directed ryanodine (Ry) -sensitive channel for Ca 2+ release from 
the jSR to the dyadic space. The jSR fluid compartment contains the Ca 2+ binding pro- 
tein calsequestin as well as the proteins triadin and junctin, which interact with the 
ryanodine receptor (RyR) and calsequestrin. This co-located configuration of the RyR 
receptor, along with the proteins calsequestrin, triadin and junctin which exist on the 
luminal side of the jSR membrane, constitutes a jSR Ca 2+ release regulating mechan- 
ism called the luminal sensor (Figure 4B; Figure 5). The protein-protein interaction 
between them plays an important role in regulating the open-state of the RyR Ca 2+ 
release channel [25]. A six-state Markovian scheme (Appendix Al, Equations 87-92) is 
used to describe the dynamics of this interaction and it is called the SR luminal Ca 2+ 
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Figure 3 Fluid compartment model. Figure 3: Representative fluid compartment model, showing 
membrane surface area separating different compartments. M FreeS ^: free SL membrane; Mj un ctionaisb 
junctional SL membrane; M yS/? : junctional SR membrane; K/\ LS r. Longitudinal SR membrane. 



sensor. Figure 4B shows a functional diagram of the luminal sensor and its output state 
is shown connected to the four-state RyR model Specifically, the sensor adjusts Ca 2+ 
dependent rate functions within the ryanodine receptor model, which affects the open 
probability P 0 of the SR Ca 2+ release channel. 

With regard to coupling, the DHP and Ry-sensitive Ca 2+ channels are assumed to be 
located on opposite sides of the small dyadic fluid space (nanodomain) as shown in 
Figure 6, and coupled functionally by a CICR mechanism. The dyadic space is assumed 
to be in fluid communication with the cell cytoplasm via a restricted diffusion region. 
In contrast, the LSR is not functionally coupled to the sarcolemma, but rather is in 
contact with the cytoplasm via the SERCA pump (as shown in Figure 3). Table 4 pro- 
vides values for parameters used to model the intracellular ion transport. 



The Dyadic Coupling Unit (DCU) 

In describing this functional unit it is necessary to provide progressively more detailed 
descriptions of the component elements of the individual dyad, particularly the geome- 
trically opposed DHP and Ry-sensitive Ca 2+ channels, as well as the geometry and 



Table 1 Surface area of various plasma membranes in the cell 



Variable 


Description 


Value 


^Ext.SL 


Surface area of external SL 


11.4 x 10 3 


Att 


Surface area of T-tubule 


5.52 x 10 3 


^TotSL 


Surface area of total SL (including external SL and T-tubule) 


16.9 x 10 3 (*) 


AjunctExt.SL 


Surface area of junctional external SL 


0.846 x 10 3 


AjunctTT 


Surface area of junctional T-tubule 


2.54 x 10 3 ^m 2 


^TotJunct 


Surface area of total junctional plasma membrane 


3.39 x 10 3 


AjunctSR 


Surface area of junctional SR 


6.99 x 10 3 


^LongSR 


Surface area of longitudinal SR 


36.8 x 10 3 


ArotSR 


Surface area of total SR 


43.8 x 10 3 



^Electrical capacitance of the cell membrane = 169 pF (using 1 juF/cm 2 ). 

Table 1: Surface area of various plasma membranes in the cell. All the above parameters are derived from Bers [16]. 
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Table 2 Parameters used to model sub-cellular morphology 



Parameter 


Definition 


Value 


References 


Ndyad 


Number of dyadic units 


10000 


[8]* 


Vmyo 


Myoplasmic volume 


5.3581 x 10~ 2 nL 


[16] + 


Vlsr 


Longitudinal SR volume 


1.1776 x 10~ 3 nL 


[16]* 


I V 

N dyad 


Total junctional SR volume 


1.104 x 10~ 4 nL 


[16]* 


Ar 


Step size in the 'r' direction 


10 nm 


Numerical solution 1 


d 


Diameter of the cylindrical cleft space in the 'r' direction 


400 nm 


[119,8,120,121,6]* 


Az 


Step size in the z' direction 


0.76 nm 


Numerical solution 1 


h 


Length of the cylindrical cleft space in the z' direction 


15.2 nm 


[119,8,120,121,6]* 


Vcleft 


Volume of a unit dyadic space 


1.91 x 10~ 9 nL 





Table 2: Parameters used to model sub-cellular morphology. Adopted (*), derived (t) or estimated (40 from the cited sources. 



buffering properties of the small dyadic space. As will be shown, we have described the 
dynamics of both types of Ca 2+ channels using Markovian state models (Figure 4) 
which include features such as Ca 2+ mediated channel inactivation, a graded CICR 
process with a "calcium gain" of approximately 6-7, and two-dimensional Ca 2+ diffu- 
sion within the dyadic space. Crank [26] discusses diffusion problems in a two-phase 
heterogeneous medium and shows that diffusion through a system of barriers (RyR 
feet structures in the dyadic cleft space) can be approximated by diffusion in the same 
region without barriers but with a reduced effective diffusion coefficient. We hence 
take this approach in modeling the Ca 2+ diffusion by solving the 2-D Laplacian 
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Figure 4 Calcium channel dynamics. Figure 4: Calcium channel dynamics, (a) Markovian model 
describing the DHP-sensitive Co 2+ channel, and (b) Markovian model of the Ry-sensitive Co 2+ channel and 
the luminal SR Ca 2+ sensor. Input from the luminal SR Ca 2+ sensor modulates the rate constants in the 
model of Ry-sensitive channel exercising the indirect bias of luminal [Co 2+ ] JS r on the RyR receptor. 
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Table 3 Parameters used to model ion transport across the sarcolemmal membrane 



Parameter 


Definition 


Value 


References 


F 


Faraday's constant 


96485 coul ■ mol 






R 


Ideal gas constant 


8314 mJ ■ mol ■ K 






T 


Absolute temperature 


290 K 


Measured 


[Co \ 0 


Extracellular Co 2+ concentration 


l .0 mM 


Measured 


[Na\ 


Extracellular No + concentration 


140.0 mM 


Measured 


[Cs ] Q 


Extracellular Cs + concentration 


3.0 mM 


Measured 




Valence of No + and Cs + ions 


1 .0 






Zca, Z Ba 


Valence ot Co and bo ions 


2.0 






Pea 


Permeability of L-Type calcium channel to Co 2+ 


6.736/ X 1 0 jL/L • S 


[30]* 




P No 


Permeability of L-Type calcium channel to No + 


8.0355 X 1 0 jL/L • 5 






Pes 


Permeability of L-Type calcium channel to Cs + 


6.2088 X 1 0 [XL ■ s 


[30]* 




KmAllo 


Dissociation constant for allosteric Co 2+ activation 


125 x 10 mm 


[85]* 




l^mCao 


Dissociation constant for extracellular Co 2+ 


1 .14 mM 


[85]* 






Dissociation constant for intracellular Co 2+ 


0.0036 mM 


[122]* 




KmNao 


Dissociation constant for extracellular No + 


87.5 mM 


[83]* 




K/nNai 


Dissociation constant for intracellular No + 


12.3 mM 


[85]* 




1/ 

v max 


K/l^virv^i i rv^ Klfi~^~ /f QvrhDnnQ r\ irront 
IVIdXIl 1 1UI 1 1 ivu /L.U cALIIdllLjc LUllcllL 


/ /O.Zj^Z LVl 


[85]* 




kmpca 


Half saturation constant for the SL Co 2+ pump 


0.5 |L/M 


[30]* 




IPMCA 


Maximum sarcolemmal Co 2+ pump current 


1.15 p/\ 


[30]* 




Kmcs 


Dissociation constant for extracellular Cs + 


1.5 x 10 3 jL/M 


[94,92,8 


8,74]* 


l^mna 


Dissociation constant for intracellular No + 


2.14 x 10 5 jL/M 


[94,92,8 


8,74]* 


iNaCs 


Maximum No + /Cs + pump current 


147.3 p/\ 


[94,92,8 


8,74]* 


^Nab 


Maximum background No + current conductance 


0.00141 nS 


[30]* 




Pa 


Mean access resistance of the tubular system 


20.0 /cD 


[123]* 





Table 3: Parameters used to model the transmembrane currents l Ca ,L , Inocw Ipmca, Inocs arid the background sodium 
current \ NQ/b . Adopted (*) or estimated (+) from the cited sources. 



equation Appendix A3 (Equations 147-150) in the DCU without explicitly accounting 
for local potential fields. The DHP-sensitive ica,L, tt channel brings in trigger Ca 2+ (0.1 
pA which is of the same order as measured by Wang et. al. [13]) causing a sparklet (a 
local increase in Ca 2+ concentration at the mouth of the channel). This trigger Ca 2+ 
causes a release from a cluster of opposing RyR channels, causing a spark. This com- 
bined release from a cluster of RyR channels causing a spark is represented as the 
release from a unitary RyR channel fo^) in our model (shown in Figure 6). The char- 
acteristics of elemental Ca 2+ release from a unitary RyR channel in our model agrees 
with data in terms of amplitude which is of the order of 3 pA (reported by Cheng 
et al. [27] and Blatter et al. [28]) and duration (full duration at half maximum 
(FDHM)) which is of the order of 50 ms (reported by Zima et al. [14]). 

The single DCU in our model represents the lumped activity of a large number of 
individual dyads (e.g. 10,000), and it is charged with the task of forming the cytosolic 
Ca 2+ transient (hence mechanical contraction) each beat of the cardiac muscle cell. In 
response to tonic application of voltage clamp pulses, the DCU strongly depends on an 
adequate supply of Ca 2+ from the SR. The measurements of Diaz et al. [29] show that, 
although trigger current may be supplied regularly by tonic voltage clamp pulses, there 
is an inherent steady-state dependence of the magnitude of Ca 2+ release on the parti- 
cular value of SR Ca 2+ content (i.e., there is a relationship between SR Ca 2+ content 
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Figure 5 Schematic of the CICR subsystem. Figure 5: Schematic of the CICR subsystem. The dyadic 
coupling unit comprises of a DHP-sensitive Co 2+ channel opposing a Ry-sensitive Co 2+ channel. The 
transmembrane proteins triadin and junction along with calsequestrin mediate the interaction between the 
luminal Co 2+ and the RyR thus regulating the release of Co 2+ flux from the jSR into the dyadic space. 



and peak [Ca 2+ ] myo ; Figure 4 Diaz et al. [29]). This plot gives us a glimpse of the 
input-output relationship of the dyadic coupling unit and indicates that SR Ca 2+ con- 
tent is an important controlling variable for the CICR process implemented by the 
DCU model. 

L-Type Ca 2+ Current 

A multiple state characterization of I Ca ,L in rat ventricular myocytes has been reported 
previously by our group [30]. The gating scheme used in this I Ca ,L model has an addi- 
tional high voltage state C6 dhpr as shown in Figure 4A which is introduced to repro- 
duce I Ca ,L tail currents. Upon voltage-dependent activation, the channel achieves the 
primary open state 02 dhpr . The degree of opening is enhanced in the presence of acti- 
vated calmodulin-dependent kinase (CaMKII act ) [31-33] which is known as Ca 1 
+ -dependent facilitation (CDF). CaMKII has been shown to tether to the Ic a ,L channel 
[34] functioning as a Ca 2+ signalling sensor for facilitation. The open probability of the 
Ica,L channel is also increased in the presence of activated calcineurin (CaN act ) [35]. 

These two effects are modeled as shown in Figure 4A, where kfj pr * s a function of 
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Figure 6 Schematic of the Dyadic nanodomain. Figure 6: Schematic of the LCC-RyR coupling. The 
dyadic cleft allows Co 2+ signaling between a DHP-sensitive Co 2+ channel opposing a cluster of Ry-sensitive 
Co 2+ channels (represented as the release channel i RyR in our model). The trigger sparklet due to Co 2+ 
entering via DHP-sensitive icaijr channel causes a spark as a result of release from a cluster of opposing 

RyR channels. The large amount of Co 2+ released subsequently diffuses into the cytosol. 

v J 

both CaMKII act and CaN act besides voltage. The gating scheme features two different 
pathways for inactivation of the open state. The pathway 02 dhpr <-> C5 dhpr accounts for 
voltage-dependent inactivation, whereas the pathway 02 dhpr <-> C4 dhpr <-> C5 dhpr with 

Qz 2+ -calmodulin (Ca 4 CaM) dependent rate constant kfl pr accounts for the fast and 

slow phases of Ca 2+ -dependent inactivation (CDI) [36]. Transitions to the closed state 
via both the inactivation pathways are suppressed in the presence of CaMKII act and 
CaN act to allow CDF. A 2-state Markovian model allows Ca 2+ mediated interaction 
between calmodulin and the l Ca ,L channel. As shown in Figure 4A, state S2 dhpn which 
denotes Ca 4 CaM bound to the IQ-motif of the Ica,L channel, modulates calmodulin 
dependent Qz 2+ -induced inactivation. This is in agreement with the findings of Nikolai 
M. Soldatov [37]. Our approach to modeling the effects of the proteins CaM, CaMKII, 
and CaN on the I Ca ,L channel was based on the premise that they are co-localized with 
the channel itself [35]. Most of the beat-to-beat modulation is produced by CaM and 
CaMKII, whereas CaN is constitutively active in the dyad [38]. Tables 5, 6 and 7 pro- 
vide values for the rate constants used to model Qz 2+ /CaM binding, CaM buffering, 
C* 2+ /CaM/CaMKII as well as Oz 2+ /CaM/CaN interactions. 
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Table 4 Parameters used to model intracellular ion transport 



Parameter 


Definition 


Value 


References 


7 PLB 


rate of PLB dp phosphorylation 


6800 s" 1 


[77]* 


7 PLB 


rate of PLB P dephosphorylation 


1 000 s" 1 


[77]* 


Kcyt,se/ra 


Maximal binding/release rate of Co 2+ from cytosol to SERCA 


6250 s" 1 


[77]* 


Kserca,sr 


Maximal binding/release rate of Co 2+ from SERCA to LSR 


6.25 s" 1 


[77]* 


SERCA tot 


Total amount of SERCA 


47.0 jjM 


[77]* 


P5R 


Phospholamban to SERCA ratio 


1 .0 


[77]* 


PKA aci 


Relative regulatory activity of PKA 


0.1 


[77] + 


Ttr 


Time const, for transfer of Ca 2+ from LSR to jSR 


7.0 X 1 0 5 


[30] 


?Na 


Time const, for transfer of No + from dyad to myoplasm 


1 .0 x 1 0 s 


[30] 


?Cs 


Time const, for transfer of Cs + from dyad to myoplasm 


6.0 X 1 0 s 


[30]* 


Pryr 


Permeability of the RyR channel 


1 "7 1/1 v ✓ 1 rx~~7 , , i -~1 
1 .714 X 1 0 jL/L • S 


r~>~7 ~>0 1 T|-t- 

L2/,28,l 3] 


D C a 


Diffusion constant for Co 2+ in the dyadic space 


1 00.0 (L/mV 1 


[1 19]* 


M, 


uciibiLy ui iiiyii diMiiiLy l.u uiiiuiiiy biicb 


?nn n /M/i 

ZUU.U (U/w 


Li i yj 


N, 


Density of low-affinity Co 2+ binding sites 


16.0 jL/M 


[119]* 


K, 


Half-saturation value of low-affinity Co 2+ binding site 


1100.0 iaM 


[119]* 


K h 


Half-saturation value of high-affinity Co 2+ binding site 


13.0 jL/M 


[119]* 


[Mg 2+ ] myo 


Intracellular Mg 2+ concentration 


634 iiM 


[124]* 


[fluo3] tot 


total concentration of indicator dye 


100.0 /jM 


Measured 


K fluo3 


association rate of Co 2+ binding to dye fluo-3 


80 /iM " 1 • s" 1 


[125]* 


kfluo3 


dissociation rate of Co 2+ binding to dye fluo-3 


90 s" 1 


[125]* 



Table 4 Parameters used to model the intracellular ion transport \ cy t,serca, ^serca.sn 'm hyR, arid Ca 2+ diffusion/buffering. 
Adopted (*), derived (t) or estimated (+) from the cited sources. 



Our previous study [30] was focused on the characterization of the I Ca ,L channel under 
conditions of low [Ca 2+ ] in the dyadic space and myoplasm. In fact, Ca 2+ - release from 
the jSR was blocked by administering a relatively high dose of ryanodine (20 ftM) in all 
experiments. Therefore, to study the additional influence of SR Ca 2+ release on Ic a ,u we 
modified the original Ic a ,L model to better characterize the process of Qz 2+ -dependent 
inactivation. In the modified I Ca ,L model, majority of the structure for the 6-state dynamic 
scheme remains the same, however changes have been made to the voltage-dependent 
inactivation rate function ( fc^f r )> an d the Ca^CaM dependent inactivation rate function 
( kfl pr )♦ The specific formulas for the Markovian state equations and the modified k^" 
and kfl pr functions used in this study are given in Appendix Al (Equations 4-50). Our 
adjustments consist of reducing the contribution from voltage-dependent inactivation pro- 
cess and strengthening the G^QzM-dependent inactivation process. The rate constants 
fodhpr anc [ fodhpr are constructed in order to provide a re-excitation window during the 
return of the clamp voltage to the resting potential. With these adjustments, the Marko- 
vian state description for I Ca>L can provide good fits to measured I Ca ,L data under both test 
conditions (presence and absence of ryanodine -sensitive Ca 2+ release) as well as produce 
tail currents during repolarization from large clamp voltages (> 40 mv). 

Calcium Buffering in the Dyadic Space 

Previous modeling work by Post et al. [39], and Post and Langer [40] considered the 
effect of Ca 2+ binding sites on the inner sarcolemmal leaflet. Following these authors, 
we included low- affinity (kj - 1.1 mM) and high-affinity (kj =1.3 ftM) Ca 2+ binding 
sites on SL wall boundary of cylindrical dyadic space. The presence of these membrane 
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bound sarcolemmal Ca + binding sites in our model has significant physiological impli- 
cations, in that it prevents local dyadic Ca 2+ concentration near the "mouth" of DHP- 
sensitive Ic a ,L channel from becoming excessively high. Allowing such a condition to 
occur can cause a reversal of the I Ca ,L current, which does not occur physically during 
normal jSR release. Addition of these SL Ca 2+ binding sites does not significantly slow 
the build-up of Ca 2+ within the dyadic space, nor is the Qz 2+ -induced Ca 2+ release 
mechanism affected [41]. 

Ca 2+ Release Channel 

The gating characteristics of the Ry-sensitive release channel are not only modulated by 
the dyadic Ca 2+ concentration at its mouth but also the jSR Ca 2+ concentration via the 
luminal sensor. Several RyR gating schemes have been deduced from isolated RyR cur- 
rents measured in lipid bilayers, including: (1) the 4-state scheme developed by Keizer 
and Levine [42] with 4 Ca 2+ binding for activation and 3 Ca 2+ binding for inactivation 
to explain the "adaptation" of the RyR observed by Gyorke and Fill [43]; (2) the 6-state 
scheme suggested by Zahradnikova and Zahradnik [44], which allowed opening the RyR 
channel upon binding of a single calcium ion; and (3) the Markov model proposed by 
Keizer and Smith [45], which can be dynamically switched among the six, five and four- 
state representations during the simulation as Ca 2+ levels vary. Stern et al. [6] demon- 
strated that none of these schemes yielded stable local control of SR release, even with 
extensive adjustment of parameters. Stern et al. [6] further reported that all schemes 
resulted in models that manifested local instability, as indicated by failure of release to 
terminate after activation, or global instability caused by spontaneous activation by rest- 
ing [Ca 2+ ] myo . Since many of the kinetic gating schemes derived from lipid bi-layer data 
fail to support stable E-C coupling in simulations, he concluded that the RyR gating pro- 
cess in situ may differ considerably from that in bi-layers. 

Our gating scheme is patterned after the release channel used in the model of Stern, 
Pizarro and Rios [46] (Figure 4B), where the channel is assumed to have four states: 
rest (closed), activated (open), inactivated (closed) and refractory (closed). The activa- 
tion gate is opened by the simultaneous, cooperative binding of two Ca 2+ ions, whereas 
inactivation depends on the binding of a single Ca 2+ ion. Four additional features have 
been built into our model: (a) the crucial role of the luminal SR Ca 2+ sensor (Figure 5) 
in assisting the inactivation of the RyR channel is modeled via the dependence of all 



Table 5 Rate constants used to model Ca/CaM binding and CaM buffering 



Rate Constant 


Value 


Rate Constant 


Value 


h CM 
^02 


4.8387* 


K 42B 


"42 


K 20 


10.0* 


K QBon 


3.5 x 10~ 4 * 


k CM 


3.4722* 


"OBoff 


1.4 x 10~ 6 * 


K 42 


500.0* 


K 2Bon 


"OBon 


K 02B 


K 02 


n 2Boff 


"OBoff 


h CM 
"20B 


fc™/100 


K 4Bon 


K 0Bon 


K 24B 


R 24 


n 4Boff 


h CM 
"OBoff 



Table 5: Rate constants used to model Ca/CaM binding and CaM buffering. Adopted (*) or estimated (t) from 
Saucerman et. al. [38]. 
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Table 6 Rate constants used to model Ca/CaM/CaMKII interactions 



Rdlfcr V_UllbLdllL 


V dlUfcr 


rvdlfc: V-UribLdlll 


VdlUc 




1.72 


h CK 
"32 


2.2 


k mPP] 


11.5 


"45 


3.35 x 10~ 3 


"21 


65.67164 


k CK 
"54 


3.4722 


k CK 
K l2 


328.3582 




2.2 x 10~ 3 


"13 


3.4722 


"64 


65.67164 


"31 


3.35 


k CK 


328.3582 


h CK 
"23 


65.67164 


h CK 
"65 


65.67164 



Table 6: Rate constants used to model Ca/CaM/CaMKII interactions. All the above parameters are adopted from 
Saucerman et. al. [38]. 



the rate constants in the 4-state RyR model, on the degree of interaction between the 
RyR and the proteins triadin and junctin (Figure 4B); (b) CaMKII act dependent 

enhancement in RyR release [47-50] is modeled via the rate functions fe[^ r , k > 

and k^ 2 r > ( c ) a stronger Ca 2+ dependent inactivation of the RyR channel (at a fixed 

depolarization level), is adopted to reflect recent observations that the inactivation of 
the RyR depends on the high local Ca 2+ concentration consequential to their own Ca 2+ 
release [51]; and (d) inactivation of the RyR at different depolarization levels is made 
dependent on local Ca 2+ concentration (i.e., [Ca 2+ ] ryr at the "mouth" of the RyR channel 
on the dyadic side) as per Wier et al. [52] and Zucchi et al. [53]. Our initial studies 

indicated that the repriming rate ( k^ 2 r ) m tne g atm g scheme of Stern et al. [46] 

was quite large, which can lead to a saturated open probability of the Ca 2+ release chan- 
nel. This occurs during the later phase of channel inactivation, where the large value of 

k^ 2 r tends to reactivate the channel, resulting in saturated Ca 2+ release. Therefore, we 

utilized a value of k 7 ^, tnat is 10% of that used by Stern et al. [46], and further assumed 

that the unitary permeation flux of the jSR release channel is proportional to jSR luminal 
Ca 2+ concentration ([C^ 2+ ] ;5 ^). The specific equations for the Ca 2+ release model are 
given in Appendix Al (Equations 73-92). 

Luminal RyR sensor 

The RyR Ca 2+ release is modulated by a multi-molecular Ca 2+ signalling complex 
which is localized to the junctional SR [54-56]. This complex consists of the ryanodine 



Table 7 Rate constants used to model Ca/CaM/CaN interactions 



Rate Constant 


Value 


Rate Constant 


Value 


h CN 
K Caon 


2.0 


h CN 
"Oon 


46.0 


n Caoff 


1.0 


h CN 
"Ooff 


537.966 




4.8387 


h CN 
"2on 


46.0 


h CN 
"20 


0.0606 


«2off 


3.2604 


"24 


3.4722 


"4on 


46.0 


"42 


0.199362 


h CN 
"4off 


1.3 x 10" 3 



Table 7: Rate constants used to model Ca/CaM/CaN interactions. All the above parameters are adopted from Saucerman 
et. al. [38]. 



Krishna et al. Theoretical Biology and Medical Modelling 2010, 7:43 
http://www.tbiomed.eom/content/7/1/43 



Page 1 5 of 66 



receptor (RyR) which functions as a Ca + conducting pore [57,58], calsequestrin (CS) 
which acts as the Ca 2+ binding protein [59,60], and the junctional SR transmembrane 
proteins triadin [61] and junctin [62]. It was known previously that the proteins triadin 
and junctin anchor calsequestrin to the ryanodine sensitive receptor [61]. More 
recently, it has been observed that the protein-protein interactions between triadin, cal- 
sequestrin and RyR modulate sarcoplasmic reticulum calcium-release in cardiac myo- 
cytes [25]. Triadin and junctin are structurally homologous proteins [63], but the 
functional differences in their roles are unclear at present. Therefore we have refrained 
from modeling these proteins separately with unique roles, and will hereafter only 
mention triadin. The RyR model proposed by Shannon et al. [64] takes a heuristic 
approach towards free SR Ca 2+ concentration dependent luminal control of the RyR 
channel. Our detailed biophysical model of the luminal sensor is strongly based on the 
recent findings of Terentyev et al. [25] which uncovers complex [Ca 2+ ]j S R) dependent, 
CS mediated mechanistic interaction of the protein triadin with the RyR channel. 

The triadin protein facilitates SR Ca 2+ release by sensitizing the RyR to activation by 
the trigger current Icu,l,tt • This is incorporated in our model by allowing the rate 
constants in the 4-state Markovian model for the RyR (Figure 4B-ii) channel to be 
functions of the activated state Al ts) which represents the degree of binding between 
triadin and RyR, in the 6-state model for the luminal sensor (Figure 4B-i).The degree 
of triadin assisted anchoring of CS to the RyR channel [61,62] is denoted by the state 
72 /5 . Triadin is also known to exist in a form bound to CS alone [25] denoted by the 
state 13 is and in its free unbound form represented as 74/ 5 in the model. CS which is 
known to exist in a Ca 2+ bound form modeled by B6 ts also modulates SR Ca 2+ release 
by influencing the open probability of the RyR channel [65-67], via interactions with 
Triadin [54,66,25]. The degree of unbound CS is denoted by B5 is in the model 

During the diastolic period the relatively large concentration of available free Ca 2+ in 
the SR results in most of the CS being bound to Ca 2+ , decreasing the degree of inter- 
action between CS and triadin. This enables a strong interaction between the available 
unbound triadin and the RyR channel increasing its propensity for trigger Ca 2+ 
induced SR release. Our model incorporates this property by facilitating movement of 
states towards Al is and B6 is in the presence of large SR Ca 2+ concentration. Following 
SR Ca 2+ release, a reduced luminal Ca 2+ concentration in the jSR causes an increase 
in the amount of free calsequestrin {BS is in the model) available to bind with triadin. 
This results in a decrease in the extent of interaction between triadin and RyR (Al is in 
the model), thus inhibiting the RyR channel and leading to robust termination of SR 
Ca 2+ release in cardiac myocytes [66]. This release termination mechanism is incorpo- 
rated in our combined RyR-Luminal sensor model (Figure 4B and Appendix Al, Equa- 
tions 73-92). 

Due to the lack of in-vivo measurements on individual state transitions, the rate con- 
stants in our novel luminal sensor model (Table 8) are chosen such that (i) the triadin 
mediated sensitization of the RyR channel (via Al ts ) provides adequate peak RyR 
release which translates into the upstroke velocity of the cytosolic Ca 2+ transient; (ii) 
the luminal sensor mediated RyR channel inactivation (via Al/ S ) causes timely release 
termination resulting in a cytosolic Ca 2+ transient duration, physiological for a rat ven- 
tricular myocyte; (iii) the rate of post-release RyR recovery results in appropriate chan- 
nel refractory characteristics. In the case of channels where there is an interaction 
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between ion flow (which is not at equilibrium) and its gating mechanism the micro- 
scopic reversibility criteria does not hold true [68,69]. The RyR channel which is solely 
modulated by [Ca 2+ ] (both on the dyadic side ([Ca 2+ ] ryr ) as well as the luminal side 
([Ca 2+ ]jSR}) experiences a strong interaction of Ca 2+ flow through itself and its gating 
mechanism described as [Ca 2+ ] ryr induced self-inhibition and the luminal sensor 
induced inactivation. Hence, the rate constants of the luminal sensor model (Figure 
4B-i) are not constrained by microscopic reversibility criteria. However, a stability con- 
straint in the form of, the sum of probabilities of all possible states corresponding to 
triadin (Al is) 12 / s , 73/ s , 74/ s ) and CS (72/ s , 13 / s , B5i s , B6i s ) being equal to one is explicitly 
imposed Appendix Al, Equations 87-88) on the model of the luminal sensor in order 
to avoid run-off. 

Ca 2+ Buffering in Myoplasm and SR 

Ca 2+ buffers play an important role in sequestering a fraction of the total Ca 2+ released 
during E-C coupling and contraction. These buffers include: (a) calmodulin (CaM), 
which is assumed to be uniformly distributed in the myoplasm and dyadic space; (b) 
troponin in the bulk myoplasm; and (c) calsequestrin (CS) in the jSR. The dyadic 
space has been shown to be accessible to calmodulin (CaM), but mostly inaccessible to 
fluorescent dyes [70,71]. Therefore, we consider the dyadic space filled with calmodu- 
lin, but not fluo-3. This provides a more direct calcium communication pathway 
between the DHP and Ry receptors. The rate constants for Ca 2+ binding to calmodulin 
were based on a model from [22], whereas those for Ca 2+ binding to troponin were 
taken from Potter and Zott [72]. Rate constants used to describe Ca 2+ binding to calse- 
questrin were based on the study of Cannell and Allen [73], whereas those for Ca 2+ 
binding to troponin-Mg complex were adopted from Lindblad et al. [74]. 

It is well recognized that fluorescent indicator dyes introduced into the cytosol also 
act as Ca 2+ buffers [22], even at submillimolar concentrations [51]. In our simulations, 
we have used 100//M fluo-3. We assume that the quantity observed experimentally as 
a "calcium concentration" signal is actually the calcium complexed with fluo-3, or 
[CaF3]. The differential equation describing the change in [CaF3] with time is given in 
Appendix A2 (Eq. 138) which follows from Shannon et al. [75]. 

Ca 2+ -Uptake 

Cytosolic Ca 2+ is pumped into the LSR (Figure 1A) by a Qz 2+ -ATPase, which lowers 
[Ca 2+ ] myo and helps to induce relaxation in cardiac muscle. The transport reaction 
involves two Ca 2+ ions and one ATP molecule [76], and it is represented by the 
description given for I C yt,serca an d I S erca,sr in Appendix Al (Equations 51-63). The model 
used for the uptake pump is adopted from Koivumaki et al. [77] and takes into 
account both the forward flux of Ca 2+ from the cytosol to the LSR lumen and the 
backward flux from the LSR to the cytosol along with the Ca 2+ buffering action of the 
SERCA protein. The phospholamban (PLB) to SERCA ratio has been fixed to 1.0, 
assuming almost equal availability of both the proteins. CaMKII act affects the SERCA 
pump via direct phosphorylation assisting in enhancement of SR Ca 2+ transport by 
increasing the pumping rate [78] and indirectly via phosphorylation of PLB [79] reliev- 
ing the inhibition caused by PLB on the SERCA pump in turn increasing the sensitivity 
of the pump for Ca 2+ uptake. These two effects are modeled by allowing the rate 
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Table 8 Parameters used in the luminal sensor model 



Rate Constant 


Value 


Rate Constant 


Value 


fe 12 


oo. I D 




^7 Q 


K 2\ 


4.1 


R 32 


2.42 


k U 
"14 


0.5 


K 35 


150.3 


K 4\ 


85.7 


K 53 


25.5 


K \2 


2.98 


K 52 


88.16 


K 43 


25.5 


K 56 


1.2 


k lS 
"34 


150.3 


k h 
"65 


401.7 



Table 8: Static rate constants used in the model for the luminal sensor. 



constants for Ca + + binding to/release from the SERCA pump as well as the rate con- 
stant for phosphorylation of PLB to be a function of CaMKII act in the cytosoL The 
activating role of CaN in modulating the SERCA pump [80,81] is accounted for in our 
model by allowing the rate constant for phosphorylation of PLB to be dependent on 
available CaN act in the cytosoL The relative regulatory role of the enzyme protein 
kinase A(PKA) is fixed at 0.1, as its modulatory effect is beyond the scope of this 
study. The current I cyt ,serca dictates the transport of Ca 2+ between the cytosol and the 
SERCA protein. Similarly, the current I serca ,sr dictates the transport of Ca 2+ between 
the SERCA protein and the LSR. The difference in these Ca 2+ currents accounts for 
the Ca 2+ buffered by the SERCA protein. The jSR is subsequently refilled by Ca 2+ dif- 
fusion from the LSR. The differential equations describing the Ca 2+ balance and parti- 
cularly the transfer between two separate SR compartments (jSR and LSR) are 
provided in Appendix Al (Eq. 72). 

Ca 2+ -Extrusion via Sarcolemmal Ca 2+ Pump 

Although the sarcolemmal Qz 2+ -pump has a high affinity for [Ca 2+ ] myo , its transport 
rate is far too slow for it to be an important factor in Ca 2+ fluxes during the cardiac 
cycle. It might, however be more important in long-term extrusion of Ca 2+ by the cell. 
Our model of the plasma membrane Ca 2+ ATPase pump current is adopted from Sun 
et al. [30]. We have used a constant value of half activation constant k mpca = 0.5/uM in 
our model on the basis of measurements by Caroni et al [82]. 

Ca 2+ -Extrusion via Na + /Ca 2+ Exchanger 

In mammalian cardiac cells, it is generally accepted that the Na + /Ca 2+ exchanger has a 
stoichiometry of 3Na + :lCa 2+ [83]. I^aCa ( a g am > the combination of InuC^tt an d InuCu, 
sl) is important in removing Ca 2+ during twitch relaxation, in competition with I up . 
A simple thermodynamic Na + /Ca 2+ exchanger current model [84] may be sufficient to 
predict the direction of Ca 2+ transport by Na + /Ca 2+ exchange and the driving force, 
however the amplitude is subject to kinetic limitations (depending on substrate con- 
centrations). A more comprehensive Na + /Ca 2+ exchanger current equation [85] is 
adopted in our model, which includes factors for allosteric Ca 2+ activation and the 
transport for Na + and Ca 2+ inside and out. The maximal flux through the exchanger 
Vmax ls estimated to ensure that the Ca 2+ ion transport (which is voltage dependent) 
via the Na + /Ca 2+ exchanger matches the influx of Ca 2+ via the I Ca ,L channel [86,87], 
maintaining whole cell Ca 2+ homeostasis. 
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Na + /Cs + Pump 

The Na + /K + pump, helps in maintaining homeostasis of the intracellular Na + ion con- 
centration. ATPase activity powers the pump, as it generates an outward Na + flux and 
an inward K + flux with a Na + to IC stoichiometry of 3 to 2 [88]. However, in our 
experimental protocol, external solution in the bath was normal Tyrode (1 mM Ca 2+ ) 
with Cs + substituted for K + in order to block the inward rectifier IC current. The 
internal solution in the pipette contained Cesium aspartate supplemented with 20 mM 
CsO, 3 mM Na 2 ATP , 3.5 mM MgCl 2 and 5 mM HEPES. 

Activation of the electrogenic sodium pump in mammalian non-myelinated nerve 
fibres [89], skeletal muscle [90] and rat brain cells [91] by Cs + is reported in the litera- 
ture. It is observed that the late effects of reducing extracellular IC concentration ([/<T + ] 0 ) 
to 0 mM in mammalian cardiac muscle can be prevented by including appropriate 
concentrations of other activator cations of the Na + IK + pump such as Cs + in the 0 mM 
[K + ] Q bathing solution [92]. Monovalent cations (including Cs + ) were also added to 
IC free bathing solution to reactivate the sodium pump in guinea-pig ventricular myo- 
cytes [93]. The effectiveness of Cs + as an external cation in activating the electrogenic 
sodium pump is known to be lesser than potassium [92]. 

We have represented IxaCs m our model by the expression for Na + IIC pump formu- 
lated by Linblad et al. [74] replacing IC ion concentrations with the Cs + ion concentra- 
tion. While ensuring whole cell Na + ion balance, the peak Na + /Cs + pump current is 
modified to be one-sixth to account for the decreased potency of the cation Cs + in 
activating the pump. The voltage-dependence of ImuCs ls adopted from the data on 
Na + llC pump from Hansen et al. [94]. 

Results 

The DCU is a fundamental element in the mechanism of CICR. The sequence of 
events resulting in CICR is triggered by the Ca 2+ entering through the I Ca ,L channel. 
The characteristics of this channel are hence examined in detail to understand its vol- 
tage and Ca 2+ dependent behavior. The ability of the l Ca ,L channel to facilitate graded 
RyR release is noted. The high gain associated with RyR release is quantified. This is 
followed by a study of the properties of the RyR channel with a particular emphasis on 
its inactivation mechanism. The Ca 2+ mediated interaction between these channels is 
investigated in detail. This is followed by an effort to understand the role of [Ca 2+ ]j S R 
in CICR. In particular, the relationship of the peak cytosolic calcium transient and the 
SR Ca 2+ content is examined. The effect of modulatory agents like caffeine and thapsi- 
gargin on SR release/uptake is studied to understand the significance of these mechan- 
isms in facilitating a normal CICR. 

L-type Ca 2+ current (l Ca L ) 

The L-type DHP-sensitive Ca 2+ channel has a key role in initiating CICR. Hence a 
thorough analysis of its activation and inactivation mechanisms is considered here. 
Though the activation of the I Ca ,L channel is solely voltage dependent once activated, 
the inactivation of the channel is influenced not only by the trans-membrane voltage 
but also the Ca 2+ concentration in the vicinity of the channel [95]. The relative con- 
tribution of the voltage-dependent and Qz 2+ -dependent inactivation pathways in the 
I Ca ,L channel model can be studied by selectively blocking each pathway in the model. 
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Figure 7 shows the model-generated waveform for I Ca>L under three conditions; (i) con- 
trol case with both voltage and calcium dependent in-activation pathways intact; (ii) 
ryanodine applied to allow only the voltage-dependent inactivation and the Ca 2+ - 
dependent inactivation via I Ca ,L self-inhibition in the absence of RyR release; and (iii) 
Ca 2+ substituted with Ba 2+ to facilitate only the voltage-dependent inactivation path- 
way and block all Ca 2+ dependent inactivation pathways (k24 dhpr set to zero). The pro- 
tocol used in case (i) was a 50 ms pulse with clamp voltages ranging from 10 mv to 40 
mv in steps of 10 mv from a holding potential of -40 mV. In case (ii) and (iii) the 
pulse duration was increased to 200 ms to replicate the data. 

Measured data obtained using the same voltage clamp pulse under all three condi- 
tions (Control, Ryan-odine application and Barium substitution) is shown for compari- 
son (data obtained from Dr. Palade's lab was pre-processed to eliminate transients 
produced by changing clamp voltage in order to obtain model fits using a non-linear 
least-squares method [18]). The degree of Ic a ,L channel opening is known to be 
enhanced in the presence of CaMKII act [33] which is known as Ca 2+ -dependent facili- 
tation (CDF). According to Bers et al. [33] the positive regulation of the channel by 
CaMKII act requires Ca 2+ influx (it is not seen when Ba 2+ is the charge carrier as in 
Figure 7D and is more strongly apparent when local Ca 2+ influx is amplified by SR 
Ca 2+ release as in Figure 7A). Increase in CaMKII act seems to play an important func- 
tion in enhancing I Ca ,L peak amplitude [96], suggesting a critical role for Qz 2+ -depen- 
dent facilitation. Our model generated results also indicate that a stronger CDF (in 
Figure 7A) in the presence of RyR Ca 2+ release causes an enhancement in Ica,L peak 
amplitude (compare the model fits to data in Figure 7A (presence of RyR release) with 
Figure 7B (no RyR release) and Figure 7C (Ca 2+ substituted by Ba 2+ )). It is important 
to note that the peak of the I Ca ,L current is reached after the RyR open probability 
reaches its maximum value owing to the faster dynamics of the RyR channel. Our 
simulations indicate that following the inward current peak, Qz 2+ -dependent inactiva- 
tion (CDI) is much faster and dominates the response for the duration of voltage 
clamp (30 ms <t <80 ms). In addition, a comparison of Figs. 7A and 7B shows that, 
Qz 2+ -dependent inactivation (CDI) caused by Ca 2+ release from Ry-sensitive channels 
is much more significant than the self-inhibition produced by Ca 2+ influx via I Ca> L 
itself. The graded behavior of the cytosolic Ca 2+ transient is shown in Figure 7C. Vol- 
tage dependent inactivation (VDI) is relatively slow compared with CDI and is best 
seen in Figure 7D where all Ca 2+ inactivation effects have been blocked by Ba 2+ substi- 
tution for Ca 2+ . Under RyR blockade (Figure 7B), the relatively slow VDI has its major 
effect during the late phase of the long voltage clamp pulse (e.g., beyond 100 ms) and 
in a time range where CDI is relatively constant. Quantitative analysis of movement of 
states (during the depolarizing pulse duration of 50ms) via different pathways in the 
six state Markovian model shows that 68.24% of the total inactivation of Icu,l,tt is via 
the calcium dependent 02 dhpr -C4 dhpr pathway, owing to the large Ca 2+ concentration 
which the channel is exposed to in the dyad. In contrast, only 19.21% of the total inac- 
tivation of Ica,L,SL is via the calcium dependent 02 dhpr -C4 dhpr pathway, because of the 
low cytoplasmic Ca 2+ concentration in the vicinity of the channel. Considering the 
total I Ca ,L current {Ic a ,L,TT + Ic a ,L,si)> around 63.34% of the I Ca ,L channel inactivation is 
via the Ca 2+ dependent pathway. 
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Figure 7 l Ca L & [Ca 2+ ] myo - comparison of model generated and experimental data in the positive 
voltage range. Figure 7: Model generated and experimentally obtained data from a rat ventricular 
myocyte clamped to voltages between 10 and 40 mv in steps of 10 mv. The holding potential used is -40 
mv. Data is obtained from Dr. Philip. T. Palade's lab. (A) l CaJL in the presence of Co 2+ release from the Ry- 
sensitive channel. The pulse duration used is 50 ms. (B) l COiL with Ryanodine administered to block RyR 
release. The pulse duration used is 200 ms. (C) Global cytosolic calcium transients [Co 2+ ] myo in the presence 
of RyR release. The pulse duration used is 50 ms (D) l Ca ,L using Barium as the charge carrier in the absence 
of RyR. The pulse duration used is 200 ms. Solid traces are the model-generated results and dotted lines 
represent the data. 

V J 



Figure 8A shows a comparison of model generated and experimentally obtained l Ca ,L 
data (obtained from Dr. Palade's lab and pre-processed to eliminate transients pro- 
duced by changing clamp voltage in order to obtain model fits using a non-linear 
least-squares method [18]) from a rat ventricular myocyte (different from the cell used 
to obtain data in Figure 7) at negative (-30 mv < v < 0 mv) clamp voltages. The corre- 
sponding graded behavior of the cytosolic Ca 2+ transient is shown in Figure 8C. Figure 
8B shows the well known [52] bell-shaped dependence of the peak [Ca 2+ ] myo on the 
clamp voltage. 

Figure 9 shows model generated normalized I Ca ,L obtained by a voltage clamp to 10 
mv from a holding potential of -40 mv using the following protocols. Case 1: Normal 
release is allowed where, following OCR, a strong Ca 2+ dependent inactivation inhibits 
the I Ca ,L current as seen in the plot numbered 1 (data obtained from Dr. Palade's lab). 
Case 2: Ryanodine applied to block RyR release; Ca 2+ entering via the lc a ,L channel 
causes Ca 2+ induced inactivation, although the magnitude of inhibition is far less than 
in the control case. Case 3: Ba 2+ substitution for Ca 2+ to completely suppress the Ca 2+ 
inactivation mechanism; the only inactivation pathway present is the slow voltage 
dependent pathway [97], which causes substantially reduced recovery compared to 
cases 1 and 2. Besides the inactivation, it can be seen from Figure 9 that the rate of 
channel activation, which is evident in the slope of the individual plots during the 
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Figure 8 l Ca ,L & [Ca 2+ ] myo - comparison of model generated and experimental data in the negative 
voltage range. Figure 8: Model generated and experimentally obtained data from a rat ventricular 
myocyte clamped to voltages between -30 and 0 mv in steps of 10 mv. The holding potential used is -40 
mv. Data is obtained from Dr. Philip. T. Palade's lab. (A) l CaJ _ invoked by negative clamp voltages in the 
range of -30 mv to 0 mv in steps of 10 mv. The pulse duration used is 100 ms. (B) Peak l Cal dependance 
on clamp voltage. (C) Global cytosolic calcium transients [Co 2+ ] myo at negative (-30 mv < v < 0 mv) clamp 
voltages. The pulse duration used is 100 ms (D) Peak cytosolic calcium transients {[Co 2+ ] myo ) dependence 
on clamp voltages. Solid traces are the model-generated results and dotted lines represent the data. 
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Figure 9 l Ca ,i_ inactivation pathways. Figure 9: l Ca ,L inactivation pathways. Voltage clamp protocol used is 
a 50 ms step pulse to 10 mv from a holding potential of -40 mv. Both model generated results and the 
data are normalized for easy comparison (1) Control case where both voltage and calcium dependent 
inactivation pathways are intact. Solid trace is the model-generated result and the dotted line represents 
the data obtained from Dr. Philip. T. Palade's lab; (2) voltage-dependent inactivation exists however Co 2+ - 
dependent inactivation is only via l Ca ,L self-inhibition without RyR release; (3) Only voltage-dependent 
inactivation pathway preserved by blocking all Co 2+ dependent pathways. 
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initial activation phase, is faster in the presence of RyR Ca 2+ release (a result of 
enhanced CDF in the presence of elevated Ca 2+ levels). This highlights the important 
role of Ca 2+ in regulating Ic a ,L channel opening and thus controlling the amount of 
extracellular Ca 2+ entering the cell 

lca f L f TT - dependent Graded SR-release 

One of the most important characteristics of the CICR mechanism is the modulation 
of the RyR Ca 2+ release based on the amount of trigger calcium delivered via the 
DHPR channel The voltage-controlled I Ca ,L channel is the link between the extracellu- 
lar excitation and the intracellular Ca 2+ release. As shown in Figure 10A, our model 
reproduces graded release. It is important to note that the onset of the Ca 2+ transient 
is also modulated based on the magnitude of the trigger calcium available to initiate 
release, as shown in Figure lOA-iii. This is due to the fact that the rate of increase in 
the open probability of the Ry-sensitive Ca 2+ release channel is controlled by the 
amount of trigger Ca 2+ present at the 'mouth' of this channel, which in-turn is graded 
by the amount of trigger Ca 2+ entering the DCU by means of the DHP-sensitive Ic a ,L, 
TT channel. This mechanism is incorporated in our model by allowing the rate con- 
stants in the 4- state Markovian model of the Ry-sensitive Ca 2+ release channel (Figure 
4) to be [Ca 2+ ] ryr dependent. Small changes in ica,L, tt cause modulation of pre-release 
[Ca 2+ ] ryn which dictates the propensity of the RyR channel for a Ca 2+ release. The 
pre-release values of the rate constants in the 4-state Markovian model of the RyR 
channel (nonlinear dependence on [Ca 2+ ] ryr as shown in Appendix Al, Equations 79- 
84) along with the diastolic [Ca 2+ ]j S R (which is kept constant in Figure 10) set the peak 
RyR open probability achieved by the RyR channel as well as the peak [Ca 2+ ] myo . The 
use-dependent adaptation of the RyR channel [13] is reflected in its non-linear 
response to trigger Ca 2+ . As shown in Figures lOB-iii, iv, the delay in the onset of cyto- 
solic Ca 2+ transient closely follows the modulation of the onset of RyR release. Though 
occurring in separate compartments, the peak of cytosolic Ca 2+ transient also tracks 
the maximum value attained by the open probability of the Ry-sensitive Ca 2+ release 
channel. Here, the clamp voltage was held constant at 10 mv to avoid the interference 
of voltage dependent Ca 2+ transport via the Na + /Ca 2+ exchanger which is co-located 
[98] with the I Ca ,L channel in the dyad and the scaling down of 1cu,l,tt corresponds to 
a fast flicker block of the channels by dihydropyridine. 

The slow rising foot that precedes the rapid upstroke (trace at position 0 in Figure 
lOB-i and B-ii) is the contribution from a single sparklet (Figure 6), which is the result 
of Ca 2+ release from a single ica,L, tt channel in our model (0.1 pA), which is of the 
same order as reported by Wang et. al. [13] bringing trigger Qz 2+ into a unitary dyad. 
This L-type Ca 2+ channel (LCC) triggers release from a cluster of RyR channels caus- 
ing a Ca 2+ spark (Figure 6), which results in the rapid increase in [Ca 2+ ] dyad that fol- 
lows the foot (shown in Figure lOB-i and B-ii). As shown in Figure 10B, the latency 
from the onset of the sparklet foot to the triggered Ca 2+ spark increases with decrease 
in the magnitude of ica,L, tt [13]. The SR Ca 2+ release flux underlying a typical Ca 2+ 
spark corresponds to approximately 2 pA [27,28,13]. From RyR single channel conduc- 
tance measurements in lipid bilayer studies [99], a Ca 2+ spark translates into a release 
from around 4 RyR channels (also reported by Blatter et. al. [28]). This single Ca 2+ 
spark corresponds to a unitary RyR release (2 pA) in our model. Our model does not 
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Figure 10 Graded CICR. Figure 10: SR calcium release graded by the trigger current Iqo.ljt ■ (A) - (i) 

Decreasing trigger current l Ca ,L,TT 00 With decreasing trigger current, RyR open probability shows a graded 
decrease in its peak as well as an increasing delay in its onset (iii) Cytosolic calcium transient follows the 
trend shown by the RyR release channel (a graded decrease in peak and increasing delay in onset) (iv) RyR 
open probability vs [Co 2+ ] ryr for increasing trigger current Ico.ljt ■ (B)- (i)Concentration profile in the dyad as 
a result of RyR release caused by a largest trigger current (see part A(i)) (ii)Concentration profile in the dyad 
as a result of RyR release caused by a smallest trigger current (see part A(i)) (iii) With decreasing trigger 
current, RyR open probability shows a graded decrease in its peak as well as an increasing delay in its 
onset (shown on an expanded time scale) (iv) Difference in concentration between the mouth of the DHP 
and Ry-sensitive Co 2+ channels. (Inset shows the cross-over when this difference becomes negative 

indicating trigger current dependent latency involved in release). 

^ ) 



attempt to reproduce the stochastic kinetics of the single channel LCC-RyR coupling; 
however it mimics accurately the average behavior of this stochastic process which is 
the net Ca 2+ release flux into the cytosol causing the Ca 2+ transient. This on/off sto- 
chastic nature of the coupling was used earlier to explain the release termination via 
stochastic attrition. However, it is now known [25] that the luminal sensor plays a fun- 
damental role in an active extinguishing mechanism [51] that effects a robust [Ca 2+ ]j S R 
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- dependent closure of the RyR channel This mechanism for inactivation of the Ry- 
sensitive Ca 2+ release channel is accounted for in our model 

Cytosolic Ca 2+ transient is also graded by the duration of the 1cu,l,tt trigger current 
controlled by the voltage clamp pulse duration. Our model also shows that triggered 
release can be prematurely stopped by rapid repolarization [100-102]. The effect of 
depolarization duration on the time course of the cytosolic Ca 2+ transient is indicated 
in Figure 11 where the duration of the voltage clamp pulse is decreased from 80 ms to 
5 ms (while keeping the clamp voltage constant at 10 mv) resulting in premature stop- 
page of release and thus a decrease in peak cytosolic Ca 2+ transient. This effect is a 
combined result of (i) the modulation of release due to pulse duration dependent 
change in the amount of trigger Ca 2+ entering the dyadic coupling unit (DCU) and (ii) 
the pulse duration dependent change in the relative role of the Na + /Ca 2+ exchanger 
co-located in the DCU. 

High Gain of Ca 2+ Release 

Besides the graded release, an extremely valuable characteristic of the CICR process is 
the high gain associ-ated with it. A small amount of Ca 2+ entering the DCU via the 
Ica,L,TT channel causes a large release of Ca 2+ from the sarcoplasmic stores via the Ry- 
sensitive Ca 2+ release channel on the jSR lumen that interacts with the DCU. This 
high gain is essential in producing physiological cytosolic Ca 2+ transients when 10000 
of these DCU's operate in tandem. Two different definitions of the gain or amplifica- 
tion factor due to CICR have been adopted in our model simulations, namely the 
ratio of: 

1. average integrated RyR flux to average integrated DHPR flux [7]; and 

2. the peak Ca 2+ transient in the presence of CICR to the peak calcium transient in 
its absence, con-tributed by the trigger calcium alone [3]. 

Thus, criterion (1) measures calcium gain observed in the dyadic space, whereas cri- 
terion (2) measures gain observed in the cytosolic space. 

With regard to the dyadic space, Figure 12A shows the Ca 2+ flux through DHPR and 
RyRs, respectively. Figure 12B shows cytosolic calcium transients under conditions of 
ryanodine blockade and with RyR Ca 2+ release. Based on our measurements, there is 
approximately a 2 ms delay from the onset of the DHPR influx to the initiation of the 
RyR release flux. 

At a clamp voltage of 10 mv, criterion (1) applied to our simulation yields an inte- 
grated flux ratio of 6.39, whereas a gain of 5.75 is calculated using criterion (2). By 
either method, a CICR amplification factor of approximately 6 is calculated, which is 
similar to that reported by Stern [3]. This result is also consistent with gain calcula- 
tions from the measured data of Fan and Palade [17] on rat ventricular cells. They esti- 
mated a gain of approximately 7 using comparisons of the rates of rise of Ca 2+ 
transients in the presence and absence of ryanodine. The model generated results are 
obtained by using a voltage clamp protocol of a 50 ms step pulse to 10 mv from a 
holding potential of -40 mv. 

The inset in Figure 12A shows the voltage dependence of CICR gain formulated 
using criterion (1). Our model shows a decline in gain as the clamp voltage is 
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Figure 11 Effect of depolarization duration on [Ca 2+ ] myo . Figure 1 1: The effect of changing the 
duration of depolarizing pulses T pd on the time course of the [Ca 2+ ] myo transient. Voltage clamp protocol 
used is a step pulse to 10 mv from a holding potential of -40 mv. Normalized [Co 2+ ] myo transient 
corresponding to depolarizing pulses of 5, 10, 20, and 80 ms duration. This model-generated result shows 
similarity to data in Figure 2, Cannell et al. [100] showing pulse-duration dependent changes in cytosolic 
Ca 2+ concentration. 



increased from -30 mv to 20 mv [52,103]. However, any further increase in clamp vol- 
tage results in a small increase in gain (Figure 1C, Altamirano et al. [104]). It is impor- 
tant to note that, with increasing clamp voltage, the decreased ability of the Na + /Ca 2+ 
exchanger (which is co-located [98] in the dyad) to extrude Ca 2+ partially compensates 
for the declining trigger current, in facilitating SR release and hence assists in increas- 
ing CICR gain. 

RyR Refractory Characteristics 

Given the fact that the CICR mechanism has a high gain associated with it, it could be 
prone to unstable behavior if not for the refractory characteristics of the Ry-sensitive 
release channel [105]. It is now known [25] that the luminal sensor plays a fundamen- 
tal role in an active extinguishing mechanism [51] effecting a [Ca 2+ ]j S R dependent 
robust closure of the RyR channel which is accounted for in our model. The RyR 
release fails to become regenerative due to the role of the luminal sensor which after a 
release occurs, forces the RyR channel into an absolute refractory state followed by a 
relative refractory state as shown in Figure 13C. When the RyR channel is in the abso- 
lute refractory period, [Ca 2+ ] ryr drops to a level much below what is caused by a spark- 
let (trigger Ca 2+ ) hence it robustly avoids re-excitation. This extremely critical 
refractory feature of the channel enabled by the RyR luminal sensor acts as a protective 
mechanism against premature Ca 2+ release in the wake of a secondary excitation that 
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Figure 12 CICR Ca 2+ gain. Figure 12: High gain associated with CICR. Voltage clamp protocol used is a 50 
ms step pulse to 10 mv from a holding potential of -40 mv. (A) Gain formulation by comparison of Ca 2+ 
flux through the DHP sensitive and RyR sensitive Co 2+ channels. The voltage dependence of gain 
measured using this method is shown in the inset. (B) Gain formulation by comparison of cytosolic calcium 
transient formed in the presence and absence of RyR Co 2+ release. 



occurs before jSR can be filled back to the control level The refractory nature of the 
channel is caused by the [Ca 2+ ]j S R dependent inhibition induced by the protein triadin 
from the luminal side of the RyR channel This restraining lock on the RyR channel 
assists in reloading the jSR. 

A dual stimulus protocol (S1-S2) was employed to study RyR refractoriness. The RyR 
channel was stimulated initially by stimulus trigger current SI, followed after an inter- 
val (T 2 in Figure 13 A) by an identical stimulus current S2. Each of the two stimuli 
directed into the dyadic coupling unit is the dyadic component {Icu,l,tt ) of trace 1 
(Jca,L evoked by a 10 mv clamp pulse) in Figure 7A with a peak amplitude of 6.657 
pA/pF and a duration (7\ in Figure 13A) of 50 ms. The stimuli SI & S2 were kept 
identical to delineate the effects of RyR refractoriness on CICR. The time interval 
between the two stimuli (T 2 ) was gradually reduced from 250 ms to 50 ms in steps of 
25 ms to observe the effects of partial recovery of the RyR channel. 

As seen in Figure 13, the decrease in interval between stimuli causes partial recovery 
of [Ca 2+ ] level in the jSR which manifests in increasing levels of unbound version of 
the protein calsequestrin, which tends to bind with triadin. As more triadin binds to 
calsequestrin, there is less left to bind to the luminal side of RyR to activate the chan- 
nel for CICR. This signaling sequence mediated by the luminal sensor via the interac- 
tion of proteins calsequestrin and triadin, allows only a partial SR Ca 2+ release based 
on the degree of SR Ca 2+ content recovery. This is an important feature that helps 
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Figure 13 RyR refractory characteristics. Figure 13: RyR refractory characteristics. (A) The stimulation 
protocol comprises of two consecutive stimuli S1 & S2, both of which are identical to the dyadic 
component {Ic3,l,tt) of trace 1 {l Ca ,L evoked by a 10 mv clamp pulse) in Figure 7 A with a peak amplitude of 
6.657 pA/pF and a duration {TJ of 50 ms. These are applied with their separation in time (T 2 ) increasing 
from 50ms to 250 ms in steps of 25 ms. (B) The RyR sensitive Co 2+ release channel shows refractory 
characteristics with an absolute refractory period of 50 ms and a relative refractory period of 250 ms. (C) 
The corresponding [Co 2+ ] JS r traces reflect partial RyR recovery. (D) Cytosolic calcium transient recovers 
completely after 250 ms. The dotted line indicates lack of RyR release when T 2 = 75 ms due to insufficient 
RyR recovery. 



restore the jSR Ca 2+ content at the end of every cycle. The RyR channel in our model 
has an absolute refractory period of 75 ms and a relative refractory period of 250 ms 
(Figure 13C). These results are consistent with the refractory period measurements by 
Sobie et al [106] on rat ventricular myocytes. When the RyR channel is in the absolute 
refractory period, [Ca 2+ ] RyR drops to a level much below what is caused by a sparklet 
(trigger Ca 2+ ) and hence robustly avoids re-excitation. However, when the RyR channel 
is in the relative refractory period, a partial release (Figure 13C) is possible, despite the 
RyR receptors affinity for release being low. 

Figures 14A and 14B show the steady state [Ca 2+ ] myo and [Ca 2+ ]j S R predicted by the 
model, with and without a functioning luminal sensor, respectively. The stimulation 
protocol used is a pulse train of amplitude (-40 mv to 10 mv), duration (50 ms) and 
frequency of 4.0 Hz which is applied for a period of 2.5 sec. The value of the luminal 
control of the RyR channel in allowing adequate SR filling can be seen in Figure 14 
where, in the presence of the luminal sensor the cell maintains a normal cytosolic 
Ca 2+ transient (Figure 14A-i). This is made possible due to sufficient SR filling (Figure 
14A-ii) as a result of RyR refractoriness. However, inactivation of the RyR channel is 
also known to depend on the high local Ca 2+ concentration consequential to its own 
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Ca 2+ release [51]. The lack of a luminal sensor forces the RyR channel to solely rely on 
Ca 2+ dependent inactivation mechanism. The resulting inadequate RyR inactivation 
depletes diastolic [Ca 2+ ]j S R level (Figure 14B-ii) causing the cytosolic Ca 2+ transient to 
diminish to new lowered steady state values (shown in Figure 14B-i). The absence of 
the luminal sensor is modeled by setting the value of the luminal sensing variable ( Var 
in the model) to a level consistent with it's normal value under 4 Hz, voltage clamp 
stimulation. This mimics the case where the luminal sensor is insensitive to changes in 
[Ca 2+ ]j S R and hence non-functional. It is important to note that the decrease (37.07%) 
in diastolic SR level (as indicated in Figure 14B-ii) is due to the lack of a luminal con- 
trol on the RyR channel resulting in a reduced rate of RyR recovery which in turn 
reflects in a compromised SR filling rate. Hence, the presence of a luminal sensor 
which monitors the SR Ca 2+ content, is key to long term Ca 2+ stability in the cell. 

CICR modulation by the jSR Ca 2+ content 

From the refractory characteristics of the RyR channel, it is evident that the RyR 
release is strongly dependent on the jSR Ca 2+ content. In fact, the SR Ca 2+ release 
through the RyR channel depends on: (a) the amount of trigger Ca 2+ entering by 
means of the Icu,l,tt channel (b) the concentration gradient of Ca 2+ between the jSR 
and the mouth of the RyR channel in the dyadic space and (c) the open probability of 
the RyR channel modulated by the interaction between the luminal sensor and the 



(ii) 



850 



1. 700 



!T 550 



U 



400 



250 



0 



Without Luminal Sensor 




500 1000 1500 2000 



500 1000 1500 2000 



(ii) 






850 




1 


700 




as 

CO 







Decreased diastolic [Ca 2+ ] jSR 



KWWVWWV 



500 



250 



2000 



1000 1500 2000 0 500 1000 1500 

Time (ms) Time (ms) 

Figure 14 Role of Luminal sensor. Figure 14: Application of a train of voltage clamp pulses (amplitude 
-40 to 10 mV; duration 50 ms) for a period of 2.5 sec. Pulse repetition rate is 4.0 Hz. Luminal sensor is 
intact in A and removed in B. Both A and B are steady state responses to voltage clamp stimulation. A-(i) 
[Co 2+ ] myo transients with normal amplitude. A-(ii) Sustained SR filling due to robust RyR inactivation. B-(i) 
Diminishing [Co 2+ ] myo transients indicate poor SR filling. B-(ii) Lowered [Co 2+ ] JS r levels indicating inadequate 
refilling of SR stores due to lack of RyR refractoriness. 
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RyR protein. The SR Ca 2+ released also inactivates the trigger current Icu,l,tt (as 
shown in Figure 7), hence causing an indirect self-inhibition of RyR release. 

Figure 15 shows a phase plot of RyR open probability (02 ryr ) versus the Ca 2+ con- 
centration at the mouth of the RyR channel in the dyadic space [Ca 2+ ] ryr constructed 
from model-generated data for a pulse of amplitude -40 mv to 10 mv and a duration 
of 50 ms. The pre-release diastolic [Ca 2+ ]j S R is 918 (aM. Following excitation by the 
trigger current during phase A, the RyR channel begins to open, allowing Ca 2+ flux 
into the dyad and elevating [Ca 2+ ] ryr at the mouth of the RyR channel. The rate of 
opening of the channel is Ca 2+ dependent; hence as [Ca 2+ ] ryr increases, RyR open 
probability increases first slowly and then ever more rapidly. With this onrush of Ca 2+ , 
[Ca 2+ ] ryr rapidly equilibrates with the [Ca 2+ ]j S R, reducing and eventually abolishing the 
concentration gradient across the channel. [Ca 2+ ] ryr soon reaches its maximum value 
(Tl in Figure 15) and begins to track the decrease in [Ca 2+ ]jSR, there being no further 
significant concentration gradient existing between [Ca 2+ ]j S R and [Ca 2+ ] ryr at the 
mouth of the RyR channel. During phase B, the rate of rise in RyR open probability 
decreases due to decreasing [Ca 2+ ] ryr and the onset of self-inhibition due to the large 
values of Ca 2+ concentration at the mouth forcing the RyR open probability to attain 
its maximum value (T2 in Figure 15). This is followed by phase C, where the RyR 
open probability begins to decrease slowly, initiating channel recovery due to decreas- 
ing overall [Ca 2+ ]d ya d levels as a result of (a) Ca 2+ fluxing out of the dyad into the 
cytosol, (b) lack of a drive from SR release due to drastically diminished Ca 2+ gradient 
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Figure 15 Ryanodine receptor open probability. Figure 15: RyR open probability modulation by 
concentration at the mouth. The pre-release diastolic [Co 2+ ] JS r is 918 jL/M. The stimulation protocol used is 
a pulse of amplitude (-40 mv to 10 mv) and duration (50 ms). Phase A - As a result of excitation by the 
trigger current, [Co 2+ ] ryr initially increases rapidly with a small accompanying rise in RyR open probability 
followed by steeper increase to peak [Co 2+ ] ryr with a large increase in RyR open probability.; Phase B - RyR 
open probability increases to its maximum value at a declining rate reflecting a falling [Ca 2+ ] ryr level and 
onset of Ca 2+ induced self-inhibition; Phase C - Decreasing Ca 2+ dyad and continued Co 2+ induced self- 
inhibition, begins the channel recovery process; Phase D - Decreasing [Co 2+ ] jSR forces a strong RyR 
inactivation via the luminal sensor hence robustly closing the channel. 
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at the mouth of the RyR channel and (c) continued Ca 2+ induced self-inhibition. 
Decreasing levels of free Ca 2+ in the jSR force a strong RyR inactivation via the lumi- 
nal sensor (beginning at T4 and continuing throughout phase D; Figure 15), which in 
turn causes a sharp decline in RyR open probability. The minimum value of [Ca 2+ ]j S R 
occurs at T5, and the channel is ultimately closed at T6 (Figure 15). The luminal sen- 
sor mediated inactivation assists in robust [Ca 2+ ]j S R recovery. 
Cytosolic peak [Ca 2+ ] dependence on SR Ca 2+ content 

Figure 16A shows three characteristic regions in the plot of the peak [Ca 2+ ] myo vs 
[Ca 2+ ]jSR. In Region I, where the jSR load is small, two things occur: (a) RyR release is 
reduced, as is the basal intracellular Ca 2+ concentration, and (b) Ca 2+ dependent inac- 
tivation of the Ica,L,TT channel is reduced due to the lowered dyadic Ca 2+ concentra- 
tion, which allows for greater Ca 2+ entry into the dyadic space via the lc a ,L channel. 
This increase in trigger current results only in a very small increase in Ca 2+ release, 
because of the inhibition of the Ry-sensitive channel caused by the luminal sensor in 




Figure 16 Peak [Ca 2+ ] myo vs [Ca 2+ ]j SR . Figure 16: (A) Relation between the peak of the cytosolic Co 2+ 
transient and the diastolic [Co 2+ ] jSR . Region I: The peak RyR open probability is very small owing to the low 
SR Ca 2+ content. This causes the peak Co 2+ transient to be very small. Region II: Increasing SR Co 2+ content 
causes an increase in the peak RyR open probability and hence in the peak of the cytosolic Co 2+ transient. 
Region III: A large rapid increase in the peak Ca 2+ myo is observed due to the large RyR Co 2+ release 
associated with the high SR Co 2+ content and saturation of RyR open probability. This model-generated 
result shows similarity to data in Figure 4, Diaz et al. [29].(B) Dependence of RyR open probability {02 RyR ) 
on diastolic Co 2+ jSR . Increasing pre-release diastolic Co 2+ jSR results in increasing peak 02 RyR . However, 
owing to the low SR Co 2+ levels the peak open probability is very small. The stimulation protocol used is a 
pulse of amplitude -40 mv to 10 mv with a duration of 50 ms. (C) Dependence of RyR open probability 
(02 RyR ) on diastolic Co 2+ jSR . Increasing diastolic Co 2+ jSR results in a linearly increasing peak 02 RyR . In the 
open probability saturation region the increase in peak cytosolic Co 2+ transient is sustained by increasing 
the area contained by each loop as seen in traces numbered 5 to 1. The progress in time occurs along the 
arrow indicated on trace 1. The stimulation protocol used is a pulse of amplitude -40 mv to 10 mv with a 
duration of 50 ms. 
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response to low SR Ca 2+ content. Thus, these factors cumulatively result in only a 
small linear increase in the peak [Ca 2+ ] myo with increasing SR Ca 2+ content in 
Region I. The behavior in Region I corresponds to traces 7-1 in Figure 16B (shows 
phase plots of 02 ryr vs [Ca 2+ ] ryr for different steady state [Ca 2+ ]j S R concentrations), 
where the peak RyR open probability shows very small gradual increase with increasing 
levels of diastolic pre-release [Ca 2+ ]j S R for a constant trigger Ca 2+ input via Ic a ,L< The 
steady state diastolic [Ca 2+ ]j S R levels being low, the maximum RyR open probability 
values attained are low. Although the peak [Ca 2+ ] ryr does not increase substantially as 
the diastolic [Ca 2+ ]j S R levels increase from 317//M to 673//M, the area enclosed by the 
[Ca 2+ ] ryr - 02 ryr loop (which indicates the amount of SR Ca 2+ released into the dyad) 
increases exponentially due to increasing peak RyR open probability combined with the 
increasing difference between rate of activation and inactivation with rate of activation 
increasing faster than the rate of channel inactivation. In Region II of Figure 16A, 
increasing SR Ca 2+ content begins to translate into a substantial increase in peak [Ca 2+ ] 
myo . This is because increasing diastolic [Ca 2+ ]j S R causes a significant increase in the 
peak RyR open probability, as shown in Figure 16C, which translates into a commensu- 
rate increase in the peak of the cytosolic Ca 2+ transient. The inhibiting role of the lumi- 
nal sensor is relieved with building [Ca 2+ ]j S R levels. The nonlinearity observed between 
the diastolic [Ca 2+ ] jSR levels 673^M (Figure 16B) and 734//M (Figure 16C) at the transi- 
tion between Regions I and II corresponds to the existence of a threshold characteristic 
for RyR release [107,108]. As the diastolic [Ca 2+ ]j S R levels increase beyond 734^M, the 
area enclosed by the [Ca 2+ ] ryr - 02 ryr loop increases substantially reaching a maximum 
value at [Ca 2+ ]j S R level of not only due to a significant increase in both peak 

[Ca 2+ ] ryr and peak 02 ryr but also due to a rapid increase in rate of RyR activation result- 
ing in a much larger difference in rate of activation and inactivation. Region III exhibits 
a large rapid increase in peak [Ca 2+ ] myo due to the large RyR Ca 2+ release associated 
with the high SR Ca 2+ content. However, it is important to note that this is not a result 
of increasing peak RyR open probability, which begins to plateau for steady state [Ca 2+ ] 
jsr values larger than 850^M (traces 6-1 in Figure 16C). The continuing rapid increase 
in peak cytosolic Ca 2+ concentration in region III is a combined result of: (a) increased 
release owing to large SR Ca 2+ content; (b) large values for saturated RyR open probabil- 
ity supported by an increase in the area contained by each loop as seen in traces 6-1 of 
Figure 16C; and (c) saturated operation of the SERCA pump which acts as a predomi- 
nant buffer in restoring the Ca 2+ concentration levels in the cytosol after RyR release 
[77]. As the diastolic [Ca 2+ ]j S R levels increase beyond 948//M, the area enclosed by the 
[Ca 2+ ] ryr - 02 ryr loop begins to decrease despite increasing peak [Ca 2+ ] ryr due to satura- 
tion of the peak RyR open probability and increasing rate of RyR channel inactivation 
due to large values of local Ca 2+ concentration {[Ca 2+ ] ryr ) assisting in faster recovery 
[51]. This model generated relationship between the peak [Ca 2+ ] myo and the pre-release 
[Ca 2+ ]j S R agrees with Figure 4 in Diaz et. al. [29]. 

As seen in traces 1-10 of Figure 16C, during the initial activation of a RyR channel 
(linear region of phase A in Figure 15), at any specific [Ca 2+ ] ryr the RyR open probabil- 
ity takes a lower value despite a higher steady state diastolic [Ca 2+ ]j S R level. This 
decreasing slope with increasing diastolic, pre-release [Ca 2+ ]j S R reflects a faster rate of 
rise in [Ca 2+ ] ryr for a larger pre-release SR content with a very small accompanying 
increase in the RyR open probability. 
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Ca 2+ Release and its Effect 

on 'Ca,L 

Caffeine 

Adachi-Akahane et al [71] investigated /c^z-induced Ca 2+ release in rat ventricular 
myocytes in the presence and absence of caffeine, which altered [Ca 2+ ] myo . Using their 
protocols and data as a guide, we modeled the effect of caffeine on jSR Ca 2+ release 
and its subsequent effect on I Ca>L . Specifically, a Boltzman relationship was adopted to 
simulate the binding of caffeine to the ryanodine receptor increasing channel conduc- 
tance, which in turn enhances jSR Ca 2+ release. The process is assumed to be dose- 
dependent, with an open probability of 0.5 at a caffeine concentration of 5 mM 
(Appendix Al, Eq. 86). 

Figures 17 A and 17B show model-generated [Ca 2+ ] myo responses in myocytes dia- 
lyzed with 2 mM Fura-2, where a simple voltage clamp pulse (-40 mV to 0 mV) is 
applied in the presence and absence of caffeine. I Ca ,L waveforms associated with these 
two protocols are shown in Figure 17C. In protocol A, conducted at normal rest levels 
of [Ca 2+ ] myo (10 nM), the depolarizing test pulse to 0 mV from the holding potential 
of -40 mV for 0.1 sec fully activates I Ca ,L > which in turn triggers a rapid but small 
amplitude Ca 2+ transient (a rise from 30 nM to 95 nM as shown in Figure 17A. Proto- 
col B starts with an application of 5 mM caffeine for a period of 1.0 sec (open state 
probability of RyR channel is 0.5), followed by the same depolarizing test pulse applied 
at 1.0 sec. The conditioning caffeine pulse induces a [Ca 2+ ] myo transient rising from a 
resting value of 30 nM to a peak value of 110 nM (Figure 17B, first peak -note the dif- 
ferent timescale). The short depolarizing pulse activates Ic a> L> which triggers a much 
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Figure 17 Effects of the application of caffeine. Figure 17: Model-simulated effects of the application of 
caffeine on \ COi ljj and [Cc? 2+ ]/ in a cell dialyzed with 2 mM Fura-2. Changes in [Co 2+ ] myo resulting from: (A) 
application of a simple depolarizing pulse from -40 mV to 0 mV, and (B) the same depolarizing pulse 
following pre-exposure to caffeine (5 mM). Panel (C) compares the \ Co ,l waveforms associated with these 
protocols. This model-generated result shows similarity to data in Figure 3, Adachi-Akahane et al. [71]. 
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smaller [Ca 2+ ] myo transient via CICR (30 nM, Figure 17B, second peak). This Ic a ,L' 
evoked Ca 2+ transient in response to the voltage pulse is smaller in amplitude than 
control (Figure 17 A), due to: (a) reduction in RyR release due to depletion in [Ca 2+ ]j S R 
([Ca 2+ ]jSR drops from 1.25 mM to 0.25 mM with caffeine application).; and (b) the 
strong competing role of the luminal RyR sensor in inactivating the RyR channel in 
response to the post-caffeine, low SR content. The /c^z-inactivation characteristics 
undergo a change as shown in Figure 17C, where trace B inactivates more slowly com- 
pared to the control I Ca ,L (trace 1). After exposure to caffeine, the overall amplitude of 
[Ca 2+ ]dh p during release does not increase as much as in the control case, since jSR 
Ca 2+ release is reduced. Consequently, Ca 2+ dependent inactivation (during release) of 
I Ca ,L a fter caffeine application is not as strong as that in the control case. However, the 
baseline level of Ca 2+ concentration at the mouth of the DHP channel {[Ca 2+ ] dhp ) is 
elevated following caffeine application, causing sustained inactivation resulting in the 
crossover of the current trace in Figure 17C. This is consistent with the observation of 
Cens et al. [41], that Ca 2+ induced inactivation is a result of calmodulin mediated sen- 
sing of the local Ca 2+ concentration. 
Thapsigargin 

Analogous to the experimental protocol of Adachi-Akahane et al. [71] dealing with 
thapsigargin blockade of the SERCA pump (their Figure 9), we apply a series of voltage 
pulses (-40 mV to 0 mV, duration 0.1 sec, frequency 0.05 Hz) to our cell model, modi- 
fied by blockade of the SERCA pump and addition of 2 mM Fura-2. Figure 18 shows 
that when Ca 2+ uptake is blocked, SR Ca 2+ content declines, as does the jSR Ca 2+ 
release with each voltage pulse (Figure 18A). Considerable time is consumed after thapsi- 
gargin application, before SR Ca 2+ content is finally exhausted (8-10 beats at 0.05 Hz; 
Figure 18A). As the magnitude of the [Ca 2+ ] myo transient decreases, the [Ca 2+ ] dhp depen- 
dent inactivation of I Ca ,L slows (Figure 18B), resulting in an enhanced peak. 

For a given clamp amplitude, a decrease in jSR Ca 2+ content decreases Ca 2+ release, 
producing decreased Ca 2+ concentration at the dyadic side of the DHP-sensitive chan- 
nel. Consequently, a smaller amount of calcium-calmodulin complex (Ca 4 CaM) is acti- 
vated, leading to a lower degree of Ca 2+ dependent inactivation. Thus, I Ca> L peak is 
enhanced, as shown in Figure 18B, increasing the area under the waveform, corre- 
sponding to an increased amount of Ca 2+ entering the cell. This inverse relationship 
between jSR Ca 2+ content and the amount of trigger Ca 2+ entering the dyadic space 
for a given clamp voltage is the essence of the feedback mechanism maintaining the 
efficacy of CICR. This mechanism is operational even in the presence of high concen- 
trations of cytoplasmic Ca 2+ -buffers, since the small dimensions of the dyadic space 
allow jSR released Ca 2+ to reach its destination (i.e., the DHP receptor) before being 
captured by the buffers [109]. 

Effect of modulation of [Ca 2+ ] Q 

Increase in [Ca 2+ ] 0 causes an increase in entry of extracellular Ca 2+ into the cell via 
the Ica,L channel. This results in an increase in SR Ca 2+ content, as shown in Figure 
19B. Aided by an increase in trigger Ca 2+ as well as a larger SR Ca 2+ content, the RyR 
release is enhanced. This increase in SR release assisted by impaired Na + /Ca 2+ 
exchange due to elevated [Ca 2+ ] 0 manifests in an increase in the peak Ca 2+ concentra- 
tion levels in the cytosol, as shown in Figure 19 A. Similarly, a decrease in extracellular 
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Figure 18 Effects of the application of Thapsigargin. Figure 18: In rat ventricular myocytes dialyzed 
with 2 mM Fura-2, thapsigargin gradually abolishes [Co 2+ ] myo transients and decreases the rate of 
inactivation of l Ca ,b This result is generated by the model, in response to a pulse train of amplitude (-40 
mV to 0 mV), duration (0.1 sec) and frequency of 0.05 Hz. (A) Time course of the effect of thapsigargin on 
A[Cc? 2+ ]j. (B) L-type Ca 2+ currents at the times indicated in (A). This model-generated result shows similarity 
to data in Figure 1 1, Adachi-Akahane et al. [71] 



Ca + levels causes a decrease in trigger Ca + available to cause CICR, which results in 
a decrease in peak [Ca 2+ ] levels in the cytosol, as shown in Figure 19C Reduced Ca 2+ 
entry into the cell via Ic a ,L> translating into decreased availability of post-release Ca 2+ 
in the cytosol, causes reduction in SR filling via the bi-directional SERCA pump. The 
corresponding drop in the SR Ca 2+ content is evident in Figure 19D. Although, our 
model-generated results show similarity to data in Figure 2, Diaz et al [110], our result 
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Figure 19 Effects of modulating [Ca 2+ ] 0 . Figure 19: The effects of changing external Co 2+ concentration 
on cytoplasmic Co 2+ and SR Co 2+ content. The stimulation protocol used is a pulse train of amplitude (-40 
mv to 0 mv), duration (100 ms) and frequency of 0.5 Hz. (A) Following sudden rise in [Co 2+ ] Q the 
instantaneous increase in [Co 2+ ] myo occurs due to more Ca entry, followed by a gradual staircase rise due 
to the the change in SR content; (B) Changes in [Co 2+ ] myo forced by a sudden decrease in [Co 2+ ] 0 ; (C) 
Changes in [Co 2+ ] jSR due to increased [Co 2+ ] OI (D) Changes in [Co 2+ ] jSR due to reduced [Co 2+ ] Q . This model- 
generated result shows similarity to data in Figure 2, Diaz et al. [110]. 



in Figure 19D does not confirm an increase in SR Ca 2+ content with a decrease in 
[Ca 2+ ] 0 as reported in Figure 2B, Diaz et al [110]. On the contrary, they agree with the 
data in Figure 5, Diaz et al 1997 which shows a decrease in free diastolic steady- state 
SR Ca 2+ content with decrease in [Ca 2+ ] 0 in cells not showing spontaneous release. In 
this study, we have attempted to isolate the effects of changes in [Ca 2+ ] 0 on intracellu- 
lar Ca 2+ levels by ensuring that there are no voltage dependent effects on Ca 2+ trans- 
port via the I Ca ,L channel, Na + /Ca 2+ exchanger or the plasma membrane Ca 2+ ATPase 
pump. The simulation protocol used in our study (reproduced from Diaz et al. [110]), 
is a pulse train of amplitude (-40 mv to 0 mv), duration (100 ms) and a frequency of 
0.5 Hz. However, it is important to note that, for any constant extracellular Ca 2+ con- 
centration, in response to a voltage clamp pulse train, increasing clamp potential (>10 
mv) results in a decrease in steady state SR Ca 2+ content as a result of the dominating 
effect of reduced Ca 2+ entry via the I Ca ,L channel despite a decrease in net Ca 2+ extru- 
sion via the Na + /Ca 2+ exchanger per cycle. 

Calcium balance under conditions of repetitive stimulation 

Our model exhibits long term Ca 2+ stability. Analogous to the experimental protocols 
of Negretti el al. [Ill], we studied the dynamic aspects of the calcium balance in the 
cell model by subjecting it to repetitive voltage clamps of different durations. For stable 
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steady-state operation, Ca + entry into the cytosol via I Ca ,L must exactly balance Ca + 
efflux. Changing the rate or pattern of stimulation can have significant effects on the 
cell's cytosolic Ca 2+ balance and subsequent contractile response [16]. To demonstrate 
Ca 2+ balance in our model, long (Figure 20A) or short (Figure 20B, C and 20D) voltage 
clamp pulses were applied at a selected repetition frequency, and the resultant changes 
in [Ca 2+ ] myo were correlated with sarcolemmal and SR Ca 2+ fluxes. Analogous to the 
protocol of Negretti et al. [Ill], pulse trains were applied from a holding potential of 
-40 mV to +10 mV, with durations of either 800 ms (long-pulse) or 100 ms (short- 
pulse) at a frequency of 0.33 Hz, for a total test duration of 2 minutes (40 pulses). In 
analyzing the Ca 2+ balance, the amount of electric charge Q transferred per pulse was 

calculated according to the equation: J T I ( t ) dt > where I(t) can be l Ca ,L > ^NaCa> Ipmca> 

I ry r> I C yt,serca> an d T is 1 cycle duration (3 ms for 0.33 Hz stimulation). 

Numerical integration was carried out with respect to steady-state levels, hence only 
activating currents and leakage currents were considered (background Ca 2+ current is 
neglected). The results of the Ca 2+ balance are presented in a manner where relative 
Ca 2+ fluxes can be easily compared. We note that in evaluating the integral of the 
combined exchanger current lNaCa> the integral is multiplied by a factor of 2 to 
account for the fact that the exchanger stoichiometry is 3Na + per Ca 2+ (the exchanger 
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Figure 20 Long term Ca 2+ stability during repetitive stimulation. Figure 20: The stimulation protocol 
used is a pulse train of amplitude (-40 mv to 10 mv) and frequency of 0.33 Hz.(A) Integrals of external and 
internal Co 2+ currents on stimulation with repetitive long pulses (800 ms); (B) Percentage of changes in 
[Co 2+ ] myo peak magnitude and cumulative SR Co 2+ loss on stimulation with repetitive short pulses (100 ms); 
(C) Integrals of external Co 2+ currents with short-pulse stimulations (100 ms); (D) Integrals of internal Co 2+ 
currents with short-pulse stimulations (100 ms). This model-generated result shows similarity to data in 
Figure 3 and Figure 5, Negretti et al. [1 1 1]. 



Krishna et al. Theoretical Biology and Medical Modelling 2010, 7:43 
http://www.tbiomed.eom/content/7/1/43 



Page 37 of 66 



transports 1 net charge per Ca + , whereas all other Ca + currents transport 2 charges 
per Ca 2+ ). 
Long-pulse protocol 

During repetitive long-pulse (800 ms) stimulation, transient Ca 2+ fluxes cross the sar- 
colemmal and SR membranes in both directions (inward (into the cytoplasm) and out- 
ward). However, [Ca 2+ ] myo returns to control levels by the end of the long-pulse 
stimulation protocol The long-pulse stimulation protocol thus serves as a means of 
studying the steady-state balance of Ca 2+ influx and efflux to and from the cytoplasm, 
under control conditions. In contrast, this Ca 2+ balance is not present in the short- 
pulse protocol, and [Ca 2+ ] myo is elevated at pulse termination. Our long-pulse simula- 
tions show that: (1) the magnitude of individual Ca 2+ transients do not change during 
the 40 pulse sequence; (2) the peak calcium release current (I ryr ) has a constant magni- 
tude; and (3) the occupancy of the calsequestrin Ca 2+ buffer in the jSR compartment is 
reduced by 15% as a result of Ca 2+ release during each cycle. Although Ca 2+ fluxes 
cross the sarcolemma and SR membranes in either direction, the integral of each of 
these currents (indicating charge transfer) over the pulse repetition interval is a con- 
stant. Figure 20A shows this well, in that, the charge transfer for each model current 
(large or small) is a constant, indicating no net loss of Ca 2+ nor consecutive-pulse Ca 2 
+ depletion (or augmentation) in the jSR compartment under the long-pulse protocol. 
In Figure 20, we refer to Ca 2+ influx and efflux through the bounding cell membrane 
(sarcolemma (SL) and T-tubule) as "external fluxes" (i.e., Ic a ,u ^NaCa> an d Ipmca)> 
whereas SR membrane Ca 2+ fluxes (I ryn I cyt , S erca) are called "internal fluxes." By con- 
vention, inward fluxes are negative (e.g. Ic a ,i) an d outward positive (e.g., I^aCa an d 
Ipmca)- Figure 20 A shows that the integrated external fluxes sum to zero. The average 
Ca 2+ charge entering the cell via Ic a ,L is 603.28 pC (pico-coulomb), and the sum of 
averaged Ca 2+ charges extruded from the cell via l^aCa an d Ipmca is 603.28 pC (602.2 
pC by luaCa an d 1.08 pC by Ipmca- Figure 20A also shows the magnitudes of the inte- 
grated internal fluxes, which sum to zero as well. In calculating inward and outward 
Ca 2+ fluxes relative to the SR lumen, we consider an inward flux (e.g., I cy t,serca) negative 
and an outward flux positive (e.g., I ryr ). The integral of I ryr is 751.8 pC, compared with 
integral of the I Ca ,L trigger current (603.28 pC). The calcium gain of CICR calculated 
in this way is low owing to the increased entry of Ca 2+ via the trigger current Ic a ,L 
for a prolonged duration of 800 ms compared to 50 ms previously. The integral of SR 
Ca 2+ uptake is 751.8 pC. These two internal component currents yield flat integral 
values during long-pulse stimulations and sum to zero. 
Short-pulse protocol 

Decreasing the length of the depolarizing pulse to 100 ms (short pulse) has a pro- 
nounced effect on the recovery of contraction (Figure 20B, C and 20D). Compared 
with the relatively flat amplitude of [Ca 2+ ] myo observed with stimulation by long pulses, 
short pulse stimulation produces a negative staircase (Figure 20B) in agreement with 
the data (Negretti et al. [Ill]), wherein the peak magnitude decays exponentially to a 
steady state level of 60% of the initial value. Figure 20B also shows that the model pre- 
dicts a similar exponential decay in jSR Ca 2+ concentration [Ca 2+ ]j S R. Both of these 
indicators show a net loss of Ca 2+ from the cell. Figure 20C and 20D show the transi- 
ent charge movements associated with the external and internal Ca 2+ currents respec- 
tively during application of the short-pulse train, which are indicative of an elevated 
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Ca 2+ load and resultant Ca 2+ imbalance when the pulse train is first applied. Balance is 
achieved, but at new lower values of [Ca 2+ ] myo and [Ca 2+ ]j S R . For convenient compari- 
son, all internal fluxes are made positive and all external fluxes negative. Figure 20C 
shows that l^aCa plays a significant role in forming the transient response seen in the 
short-pulse protocol in extrusion of excess Ca 2+ load on the cell. The decreasing 
amplitude of the Ca 2+ transient causes a decreased LSR uptake current {I cy t,serca)> 
which in turn causes a decrease in [Ca 2+ ]j S R resulting in decreased release {I ryr ). With 
successive pulses, I cy t,serca and l ryr decline further, mirroring the exponential decline in 
luaCa toward balance. At steady state, the external and internal fluxes reach the new 
equilibrium, where the sum of the integrals of external currents and the sum of the 
integrals of internal currents are zero. 
Long-term calcium stability at higher pacing rates 

Most VC experiments on cardiac cells are usually conducted at low stimulation rates. 
Our simulations of typical VC experiments at these low rates all exhibit long term cal- 
cium stability (Figures 17, 18, 19 and 20). However, our model can also exhibit sus- 
tained calcium balance at higher pacing rates. Figure 14A shows steady state cytosolic 
Ca 2+ transients in response to a repetitive 4 Hz (which is more physiological for a rat 
ventricular myocyte) voltage clamp stimulation lasting for 2.5 seconds. The sustained 
peak and basal Ca 2+ levels of [Ca 2+ ] myo over this prolonged time frame are indicative 
of long tern Ca 2+ stability in the model. The corresponding [Ca 2+ ]j S R profile over 2.5 
seconds is also shown in Figure 14A, which indicates sustained SR filling owing to the 
luminal sensor mediated control of the RyR channel. The stimulation protocol used is 
a pulse train of amplitude (-40 mv to 10 mv), duration (50 ms) and frequency of 4.0 
Hz. Tables 9 and 10 provide values for the initial conditions used in our model. 

Secondary [Ca 2+ ] myo transients induced by "tail currents" 

Our voltage clamp simulations have dealt with clamp voltages in the range -30 mV < 
V < 40 mV. Some experimental studies [100,2,112] however have employed even larger 
clamp voltages (40 mV < V < 60 mV) to explore the high voltage behavior of the 
DHP-sensitive calcium channel. Such large clamp voltages elicit a brief Ca 2+ influx 
called a "tail current", which has been shown to trigger RyR release and hence cause 
contraction during repolarization [100,2]. The secondary Ca 2+ transient induced by the 
"tail current" is a critical argument in favor of CICR as these "tail currents" are not 
observed in skeletal muscle where the membrane potential directly controls SR Ca 2+ 
release [112]. Our model reproduces this secondary Ca 2+ transient observed during the 
return to resting potential from a large clamp voltage (> 40 mV). Figure 21 A shows a 
cartoon depicting the voltage clamp stimulation protocol used in our study where a 
pulse of amplitude (-40 mV to +50 mV) is employed with the pulse duration (T p ; Fig- 
ure 21A) increasing from 50 ms (trace 1) to 200 ms (trace 7) in steps of 25 ms. The 
clamp potential transition time (T t ) is fixed at 1 ms. The peak of the "tail current" at 
the end of the pulse is 25 times larger than the peak of the I Ca ,L current observed at 
the beginning of the pulse. This is in agreement with model generated VC data 
reported by Geenstein et al.([8]; Figure 8). Figure 21B shows that, with an increase in 
clamp pulse duration, there is a corresponding increase in the peak of the trail current 
(see traces 1-7) until the effect ultimately saturates. Corresponding to this increase in 
peak Ica,L ta il current, there is an increase in the open probability of the RyR Ca 2+ 
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Table 9 Initial Conditions 



State variable 


Definition 


Initial Value 


[Co ] myo 


Co 2+ concentration in the myoplasm 


8.1027 x 10 mM 


[Co \jSR 


Co 2+ concentration in jSR 


1 ■">/■ — 7 ~1 A A 

1 .26/7 mM 


[Co \ LSR 


Co 2+ concentration in LSR 


1 .3346 mM 


[No + ] myo 


No + concentration in myoplasm 


1 r ~7 a r a a 

1 6.746 mM 


[Na + ] dyad 


No + concentration in dyadic space 


16.321 mM 


l>-5 \myo 


Cs + concentration in myoplasm 


140.2154 mM 




Cs + concentration in dyadic space 


140.2157 mM 


o c 


Fractional occupancy of Calmodulin by Co 2+ 


0.033091 


Ore 


Fractional occupancy of Troponin by Co 2+ 


0.016049 


Ctmgc 


Fractional occupancy of Troponin by Mg 2+ and Co 2+ 


0.321 764 


Ctmgmg 


Fractional occupancy of Troponin by Mg 2+ 


0.598385 


[CoF 3] 


Concentration of Fluo3-Q? 2+ complex 


21.88721 jjM 


C1 ry r 


Closed (resting) state of RyR channel 


0.9990953 


02 ryr 


Open (activated) state of RyR channel 


1 .03668 X 1 0 


C3 ryr 


Inactivated state of RyR channel 


r\ 0011 1 w 1 r\~ 1 3 

9.3871 1x10 


01 dhpr 


Closed (resting) state of DHPR sensitive Co 2+ channel 


0.1 673614 


02 dhpr 


Open (activated) state of DHPR sensitive Co 2+ channel 


1 A C\C\ 1 "7") v / 1 (~\~^> 

1 .4991 73 X 10 


03 d hpr 


Open (activated) state of DHPR sensitive Co 2+ channel 


3.300291 x 10 


C^dhpr 


Closed (resting) state of DHPR sensitive Co 2+ channel 


7.478058 x 10" 8 


C6 dhpr 


Closed (resting) state of DHPR sensitive Co 2+ channel 


7.478058 X 10 


A] is 


Fraction of Tr/J bound to RyR 


0.6709816 


His 


Fraction of Tr/J bound to RyR and Calsequestrin 


0.1 55258 


Bis 


Fraction of Tr/J bound only to Calsequestrin 


0.0623695 


B6, s 


Fraction Calsequestrin bound to Co 2+ 


0.6 1 9346 


[Co ] serca 


Ca bound to the serca protein 


36.0637 X 1 0 jdM 


52-dhpr 


Fraction of IQ Motif bound to Ca 4 CaM 


6.816729 X 10 


PLB dp 


Fraction of unphosphorylated Phospholamban 


7.684160 X 10 


Co 2 CoM 


2 Co 2+ ions bound to C-terminus of CaM 


~> a urn ~>n 

34.56529 


C04C0M 


4 Co 2+ ions bound to C & N terminus of CaM 


8.635052 X 10 




Ri iffprprl CxKA 


/ JOjOjO a I u 


C0 2 C0MB 


Buffered Co 2 CoM 


2.035086 


Co 4 CoMB 


Buffered Co 4 CoM 


1.288455 x 10" 6 


Co 4 CoN 


CaN bound to 4 Co 2+ ions 


2.606246 x W A iiM 


C0MC0N 


CaN bound to CaM 


4.348535 x 10" Vw 


C02C0MC0N 


CaN bound to 2 Co 2+ ions and CaM 


1.419613 x 10" V/W 


Co 4 CoMCoN 


CaN bound to 4 Co 2+ ions and CaM 


3.47341 2\iM 


[Co 2+ ] dyad 


Co 2+ concentration in the dyadic space 


9.012 x 10" 5 mM 



Table 9: Initial conditions used for the state vector in order to solve the linearized system of ordinary differential 
equations. 



channel, with RyR Ca + release indicated by the corresponding decrease in [Ca + ]j S r 
(panel C) and increase in [Ca 2+ ] myo (panel D). We note from Figure 21C that for a 
shorter 50 ms VC pulse, there is no RyR Ca 2+ release although an I Ca ,L tail current is 
produced. This tail current produced with the 50 ms pulse slightly augments the 
[Ca 2+ ] myo transient just after its peak (trace 1, panel E), but there is no corresponding 
decrease in the [Ca 2+ ]j S R transient (trace 1, panel D) indicating no RyR Ca 2+ release. 
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Table 10 Initial Conditions 



State variable 


Definition 


Initial Value 


Pi 


Fraction of inactive dephosphorylated CaMKII in Co 2 CoM bound state 


5.527608 x 10~ 1 


Pi 


Fraction of active dephosphorylated CaMKII in Co 4 CoM bound state 


3.661260 x 10~ 1 


Pe 


Fraction of active "T/ir^-autophosphorylated but CaM autonomous CaMKII 


1.314410 x 10~ 3 


P 5 


Fraction of active 7/i/^ 87 -autophosphorylated but Ca2CaM bound CaMKII 


6.277911 x 10~ 7 


PA 


Fraction of active 7/i/^ 87 -autophosphorylated but Ca4CaM trapped CaMKII 


9.121920 x 10~ 8 



Table 10: Initial conditions used for the state vector in order to solve the linearized system of ordinary differential 
equations. 



At the end of this short pulse (50 ms), the RyR channel is in fact in its absolute refrac- 
tory period. The small secondary increase in [Ca 2+ ] myo transient as seen in data corre- 
sponding to trace 4 of Figure 7C is a result of the large Ca 2+ influx via the I Ca ,L "tail 
current". During this phase, the Na + /Ca 2+ exchanger is biased to extrude Ca 2+ out of 
the dyad, and hence cannot be responsible for this secondary rise in the cytosolic Ca 2+ 
transient. As pulse duration is increased beyond 100 ms the "tail current" causes a 
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Figure 21 Secondary [Ca 2+ ] myo transients induced by "tail currents". Figure 21: The stimulation 
protocol used is a pulse of amplitude (-40 mv to +50 mv) with pulse duration T p increasing from 50 ms 
(trace 1) to 200 ms (trace 7) in steps of 25 ms. (A) A cartoon depicting the voltage clamp stimulus protocol 
where, pulse duration (Tp) is varied to obtain traces 1-7 and the clamp potential transition time (7 f ) fixed at 
1 ms. (B) l COiL tail currents elicited during repolarization from +50 mv to -40 mv. (C) Open probability of the 
RyR channel indicating gradual recovery from the refractory period (D) [Co 2+ ] JS r traces show increase in RyR 
release with increasing pulse duration as a result of adequate RyR channel recovery. (E) Cytosolic calcium 
transients indicate secondary RyR release caused by tail currents. 
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gradual increase in RyR open probability (traces 3-7 in Figure 21) indicating progres- 
sive recovery from the refractory period. Minimum value of [Ca 2+ ]j S R in traces 3-7 of 
Figure 21D show increase in SR Ca 2+ release with increasing pulse duration as a result 
of adequate RyR channel recovery and increasing I Ca ,L "tail current". Cytosolic Ca 2+ 
transients shown in Figure 2 IE indicate "tail current" induced secondary RyR release 
for pulse duration > 100 ms. It has also been previously reported [113] that recruit- 
ment of an additional population of previously 'silent' Qz 2+ channels could cause facili- 
tation of tail currents at increasingly large clamp voltages with a time-dependence 
associated with the recruitment process. Our model not only accounts for the contri- 
bution of the already open channels (in states 02 dhpr tk03 dhpr ) to the tail current dur- 
ing repolarization but also allows for the voltage- and time-dependent recruitment of 
'silent' Ca 2+ channels modeled using a high voltage state C6 dhpr . This study proves to 
be an adjunct to the study presented in Figure 13 in establishing the refractory nature 
of the RyR channel. 

Cytosolic Buffering 

In our model, cytosolic buffering is attributed to several factors including: (a) calmodu- 
lin (CaM); (b) the C^ 2+ -specific (Tc) troponin binding site; (c) the Ca 2+ - Mg 2 * compe- 
titive troponin binding site; and (d) the fluorescent indicator dye Fluo3 used to detect 
changes in [Ca 2+ ] myo . The effects of the component buffers in helping to maintain a 
low [Ca 2+ ] myo are shown in Figure 22 in terms of occupancy functions, such as O c 




Time (ms) 

Figure 22 Occupancy of Ca 2+ buffers. Figure 22: Contributions of fractional occupancy of Co 2+ buffers in 
SR and cytoplasm. Transfer from 0 C and 0 tc to 0 tmgc occurs as a result of the pulse stimulation. 0 tmgmg is 
shown as a dotted line as it does not reflect Co 2+ binding. The corresponding [Co 2+ ] myo is overlayed (note 
the different axis) 

^ ) 
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(fractional occupancy of calmodulin by Ca + ), O tc (fractional occupancy of troponin-Ca 
sites by Ca 2+ ), O tmgc (fractional ocupancy of troponin-Mg sites by Ca 2+ ), O tmgmg (frac- 
tional occupancy of troponin-Mg sites by Mg 2 *) in the cytosol, and O ca i se (fractional 
occupancy of calsequestrin by Ca 2+ ) in the jSR Ca 2+ release compartment. Figure 22 
shows that when the Ry-sensitive Ca 2+ release channel is triggered, the jSR releases 
Ca 2+ and the occupancy O ca i se in the jSR, declines from 68% to 37%. Ca 2+ release from 
the jSR induces fast Ca 2+ binding by calmodulin and troponin in the cytosol, as repre- 
sented by the increases of O c (0.29) and O tc (0.16) in Figure 22. Interactions between 
Ca 2+ and the troponin-Mg sites result in an increase in occupancy of these sites by 
Ca 2+ (O tmgc ), and a decrease in occupancy by Mg 2 * {O tmgmg ). 

Discussion 

It is well-established that mammalian cardiac excitation-contraction coupling is 
mediated by calcium-induced calcium release (CICR). We have developed a compre- 
hensive mechanistic model of CICR under voltage clamp conditions in the rat ventri- 
cular myocyte, which includes electrical equivalent circuit models for both the free 
sarcolemma and that portion involving junctional transmission, as well as, fluid com- 
partment models for several fluid media within the cell (dyadic cleft space, longitudinal 
sarcoplasmic reticulum (LSR or Ca 2+ uptake compartment), junctional sarcoplasmic 
reticulum (jSR or release compartment), and the cytosolic fluid compartment). An 
external bathing medium completes our fluid compartment description of the cell Fig- 
ure 1. The multiple component model is referred to as the "whole-cell model" (Figures 
1 and 2). We have probed the mechanisms regulating Ca 2+ in the cell. In particular, 
we have focused on the dyadic mechanisms effecting CICR. 

The dyadic controller is a finely tuned coupling device consisting of two opposed 
Ca 2+ channels separated by a small dyadic space (Figure IB): the sarcolemmal DHP- 
sensitive "trigger" channel and the Ry-sensitive jSR Ca 2+ release channel. The trigger 
channel is voltage-activated and is driven in our simulations by a voltage clamp 
pulse which opens the channels, admitting Ca 2+ influx (trigger current) to the dyadic 
space. After diffusion within the dyadic space, trigger Ca 2+ affects the Ca 2+ depen- 
dent open probability of the jSR Ca 2+ release channel over a small range of Ca 2+ 
concentrations. Clamp voltage magnitude strongly affects the I Ca ,L current and hence 
the amount of Ca 2+ delivered to the dyadic space, but the range over which the 
Ca 2+ concentration at the mouth of the Ca 2+ release channel changes (in the 
absence of RyR release) is quite small (23.42^M to 8.82^M for a change in VC from 
10 mv to 40 mv). The RyR-mediated Ca 2+ release is activated by trigger current, but 
the release itself is affected by the concentration gradient between the jSR and the 
dyadic space, and the temporal and refractory properties of the ryanodine receptor 
(Figure 13). The four-state RyR model is informed regarding supply Ca 2+ content by 
the jSR luminal sensor, a novel feature in our model which characterizes the impor- 
tant protein-protein interactions between calsequesterin, triadin/junctin and the RyR 
receptor. Triadin/junctin strongly regulates the sensitivity of RyR to trigger Ca 2+ . 
Hence the luminal sensor, which is a key element responsible for robust post-release 
RyR inactivation and refractoriness of the Ry-sensitive Ca 2+ release channel, is criti- 
cal in providing realistic fits to cytosolic Ca 2+ transients and an adequate refilling 
time for the SR Ca 2+ stores. Ca 2+ release is thus graded with Ca 2+ concentration at 
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the mouth of Ca 2+ channel in a very sensitive manner with a gain of approximately 
7, as is shown in Figure 12. 

The sarcolemmal portion of the dyadic membrane defining the dyadic space contains 
voltage-sensitive Ca 2+ channels and deals with changes in the external environment of 
the ventricular cell (e.g., membrane response to changes in transmembrane potential 
and chemical signaling agents). It integrates these various stimuli and delivers a trigger 
current to the small dyadic space. In contrast, the jSR membrane lining the opposite 
boundary of the dyadic space is concerned with the adequacy of Ca 2+ release in CICR. 
It contains Ry-sensitive Ca 2+ channels that require a Ca 2+ concentration gradient 
directed across the channel and into the dyadic space for their operation. Thus, jSR 
Ca 2+ concentration must be maintained within an acceptable range so that calcium is 
always available for ready release. The relationship between jSR Ca 2+ content and peak 
[Ca 2+ ] myo is shown in Figure 16 A. The LSR compartment connects the jSR compart- 
ment with the cytosolic compartment. It feeds make-up calcium to the jSR, by using a 
SERCA pump to actively draw in Ca 2+ from the myoplasm. Pumping rate is controlled 
by the cytosolic Ca 2+ concentration [Ca 2+ ] myo as well as the LSR Ca 2+ concentration 
[Ca 2+ ] LS R as part of a mechanism for replenishing and maintaining jSR Ca 2+ stores. 
Our model also incorporates the Ca 2+ induced CaM mediated effects of CaMKII and 
CaN on targets such as the DHP-sensitive Ic a ,L channel, Ry-sensitive Ca 2+ release 
channel, as well as the SR Ca 2+ ATPase pump. Provision for the effects of phospho- 
lamban on SERCA has also been included. 

The voltage-gated and chemically gated channels of the dyad are tightly coupled by 
feedback mechanisms that involve Ca 2+ signaling. Although physically separated, the 
voltage-sensitive Ca 2+ channel (Figure 9) as well as the Ry-sensitive release channel 
(Phase C in Figure 15) are inhibited by increased Ca 2+ levels in the dyadic space [53]. 
In earlier models which lacked the luminal sensor, RyR self-inhibition by dyadic Ca 2+ 
was the only mechanism besides stochastic attrition (both of which are inadequate in 
effecting a robust RyR channel closure) that was given the role of initiating RyR recov- 
ery. This has been tested in our model by artificially clamping all the variables (includ- 
ing concentration at the mouth of the RyR channel) to the levels reached at T2 in 
Figure 15 despite which the RyR open probability, begins to decline in phase C of Fig- 
ure 15 showing the self-inhibitory role of high dyadic Ca 2+ concentration. Our model 
predicts (data not shown) that this Ca 2+ signaling continues even in the presence of 
high concentrations of Ca 2+ buffering agents in the cytosol (in agreement with the 
data of Adachi-Akahane et al. [71] and Diaz et al. [110]). Tight coupling between the 
DHP and Ry-sensitive Ca 2+ channels within the dyadic space thus preserves the 
mechanism of CICR under these extreme conditions. 

Regulation of cytosolic Ca 2+ concentration [Ca 2+ ] myo is evident at the cellular level. It 
is affected strongly by the voltage and [C^ 2+ ] mjo -dependent properties of the sarcolem- 
mal currents l^aCa and Ic a ,u as well as the [C^ 2+ ] mjo -dependent properties of the 
Ipmca and SERCA pumps. We demonstrate a model-generated whole-cell Ca 2+ bal- 
ance, which shows the importance of the Na + /Ca + exchanger in extruding the Ca 2+ 
that has entered the cell under normal activity, and also any excess that might occur 
when cytosolic Ca 2+ levels rise. The variety of experiments emulated in this study 
demonstrates quantification of Ca 2+ balances for all external and internal Ca 2+ fluxes 
and shows that the model has long-term stability in regulating cytosolic Ca 2+ , as 
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shown in the 120-sec duration experiments of Negretti et al. (Figure 20) at a 
pulse repetition rate of 0.33 Hz., and the faster paced stimulation at 4 Hz shown in 
Figure 14A. 

We have examined the dyadic component of the Ic a ,L with regard to its Ca 2+ 
-dependent inactivation as a function of jSR Ca 2+ content (Figure 17C). Trigger cur- 
rent is voltage-dependent and therefore parameterized by clamp pulse amplitude (Fig- 
ure 7A). In addition, its inactivation is Ca 2+ dependent, especially during Ca 2+ release. 
At a constant depolarizing voltage pulse, the pre-release jSR Ca 2+ content dictates the 
peak of the cytosolic Ca 2+ transient (Figure 16A), by controlling peak Ca 2+ concentra- 
tion in individual dyads. As jSR Ca 2+ content increases, so does the peak dyadic Ca 2+ 
concentration, and the amount of sarcolemmal Ca 2+ influx declines due to greater 
Ca 2+ dependent inactivation. This autoregulatory feedback mechanism helps to estab- 
lish a stable operating point for jSR Ca 2+ content. Transient increases in jSR Ca 2+ con- 
tent bring about increases in Ca 2+ release, but reflexly, decrease sarcolemmal Ca 2+ 
influx via Ic a ,L- Coupled with other dyadic and extra-dyadic mechanisms {luaCa) that 
decrease [Ca 2+ ] myo and hence Qz 2+ -uptake via the SERCA pump, jSR Ca 2+ content 
decreases. For decreases in jSR Ca 2+ content, the opposite occurs, so that regardless of 
the sign of the perturbation in jSR Ca 2+ content, it tends to stay constant in the 
steady-state. This is an important feature of the dyadic mechanism that preserves the 
integrity of CICR. The gain of this feedback system about the operating point is visua- 
lized as the slope of the peak Ca 2+ transient vs SR Ca 2+ content characteristic (Figure 
16 A) for a given voltage pulse level. Observing this figure, we note that there is a linear 
operating range beyond which the system gain increases dramatically. Our model 
results also show that a decrease of jSR Ca 2+ content (simulated either by a decrease 
of Ca 2+ uptake into the LSR by thapsigargin or an increase of Ca 2+ leak out of the jSR 
by caffeine) decreases systolic [Ca 2+ ] myo , and hence the model might serve as a useful 
adjunct in a study of heart failure, where decreased contractility as a result of dimin- 
ished Ca 2+ transients are commonly observed [114]. 

Model Limitations 

(a) This model of a rat ventricular myocyte is limited to Ca 2+ related channel, 
exchanger and pumps {Ic a ,u ^NaCa> Ipmca and SERCA pump), while lacking exclu- 
sive Na + or IC related channels and transporters and is based on data at positive 
potentials in the range 10 mV < V < 40 mV. It is aimed at mimicking voltage 
clamp conditions where channels other than calcium are blocked, and it cannot be 
used to study experiments involving the generation of action potentials. However, 
its focus on the Ca 2+ dynamics allows one to comprehend more clearly the impor- 
tant role of Ca 2+ signalling pathways and feedback control systems in maintaining 
whole cell homeostasis over a prolonged period of time. 

(b) A single dyadic space in our model has one representative, lumped DHP-sensitive 
and Ry-sensitive Ca 2+ channel, on opposing sarcolemmal and SR surfaces respec- 
tively. This simplified configuration does not conform to detailed structural informa- 
tion regarding the geometrical relationships between DHP and RyR-sensitive Ca 2+ 
channels [115,55,116] and thus cannot be used to draw conclusions about this part 
of the EC coupling process. However, our lumped abstraction forms a functional 
model of the dyadic coupling unit that can produce accurate predictions of cytosolic 
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Ca + transients. The effectiveness of this model is further demonstrated by its ability 
to accurately characterize the interaction between the DHP and Ry-sensitive Ca 2+ 
channels, including pulse duration dependent termination of release (Figure 11), Ca 2 
+ dependent inactivation of I Ca ,L (Figures 7, 9), as well as the wide variety of whole- 
cell voltage clamp protocols (Figures 17, 18, 19 and 20). 

(c) Although our model provides secondary Ca 2+ tail transients elicited by I Ca ,L "tail 
currents" (Figure 21), this aspect of the model has not been verified extensively due 
to paucity of measured VC data showing tail transients over the high voltage range 
(40 < V < 60) in rat ventricular myocytes. Our lc a ,L tail currents are in general 
agreement with model generated data at a clamp voltage of 50 mV reported by 
Greenstein et al. ([8]; Figure 8). The restitution time for the RyR channel in rat 
ventricular myocytes is believed to be at least 25 ms [117] and as fast as 150 ms 
[118] indicating the wide range of values reported in the literature. We demon- 
strate model-generated RyR refractoriness by changing the duration of the simu- 
lated VC pulse and obtain RyR recovery characteristics that are consistent with 
measured data for Ca 2+ spark restitution in rat ventricular myocytes reported by 
Sobie et al. [106]. However, our results on this important phenomena, particularly 
the onset of tail transients with increasing inter-stimuli interval are preliminary and 
further modeling investigations would benefit considerably from the availability of 
additional measured data. 



Conclusion 

We have developed a mathematical model of Ca 2+ dynamics under voltage clamp 
conditions in the rat ventricular myocyte, which is based solidly on experimental 
data and includes the most extensive description available of a novel feature, namely 
the luminal Ca 2+ sensor in the junctional SR which models the protein-protein inter- 
action between triadin/junctin, calsequestrin and the RyR channel. The luminal sen- 
sor imparts the much needed refractoriness to the Ry-sensitive Ca 2+ release channel. 
This element is critical in providing realistic fits to cytosolic Ca 2+ transients and an 
adequate refilling time for the SR Ca 2+ stores. Our voltage-clamp simulations 
demonstrate graded Ca 2+ transients with sufficient gain, as well as quantification of 
Ca 2+ balances for all external and internal Ca 2+ fluxes. Our model of the dyadic cou- 
pling unit (DCU) provides mechanistic explanations of the major input-output rela- 
tionship for CICR (Figure 16), as well as its modulation by trigger current (clamp 
voltage). The variety of experiments emulated in this study demonstrates that the 
model has long-term stability in regulating cytosolic Ca 2+ , as shown in the 120-sec 
duration experiments of Negretti et al. (Figure 20) at a pulse repetition rate of 0.33 
Hz., and the faster (physiological) paced stimulation at 4 Hz shown in Figure 14A. It 
also provides biophysically based insights into the molecular mechanisms underlying 
whole-cell responses to the wide variety of testing approaches used in voltage clamp 
studies of myocytes that have appeared in the literature over the past two decades 
(Figures 17, 18, 19 and 20). Thus, the model serves as a platform for the predictive 
modeling of VC investigations in a number of areas. These include new hypotheses 
with regards to the under-expression of triadin/junctin resulting in a malfunctioning 
luminal sensor, which could affect long-term calcium stability of the cell (Figure 
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14B), and/or changes in the refractoriness of the RyR Ca 2+ channel (Figures 13 and 
21) affecting the integrity of CICR under a variety of conditions. These are funda- 
mental issues that would benefit from a better mechanistic understanding of 
deranged calcium signalling in the rat ventricluar myocyte. This study is aimed at 
providing an initial step towards this goal. 

Appendix 

Below is the complete set of equations used in the model. 

Al - Equations for currents in the model 
L-Type Ca 2+ current 
Ca 2+ current through the DHP-sensitive I Ca ,L channel 



F 2 V 



^Ca,L ~ R CaA 02 dhpr + 03 dhpr) P Ca Z Ca X 

2FV 

[Ca 2+ \ dhp e^ -Ml.0[Ca 2+ \ 0 

2FV 

-I 

Na + current through the DHP-sensitive I Ca ,L channel 

( v-21.73 ) 
21.23841 J 



1.O56X0V 



l Ma v-21.73 

1.056 xe 21 23841 + [Ca 2 0 + ] 

' FV 

0.75[Na + ] dhp eRT -0.75[Na + ] 0 

FV_ 
o RT _1 



F 2 V 



(1) 



(2) 



Cs + current through the DHP-sensitive I Ca ,L channel 



1.056xe 



v-21.73 
21.23841 



l Cs 



( v-21.73 ) 
I 21.23841 J + [Ca 



" Pcs^Cs 



F 2 V 
RT 



1.056xe 

( FV 

0.75[Cs + ] dhp e^ -0. 75 [Cs + ] 0 

FV 

e k Cs RT _j 



(3) 



where k Cs = 0.5 and [Ca 2+ ] dhp , [Na + ] dhp , [Cs + ] dhp are Concentrations at the mouth of 
the DHP-sensitive Ca 2+ channel. 

The corresponding unitary currents ic a >L> iNa> ic s are obtained by dividing the above 
net channel currents by the number of dyadic units, N dyad , 

Gating scheme for the 2-state Markovian model used to allow Ca 2+ mediated inter- 
action of the DHP-sensitive I CafL channel and calmodilin: 



SI 



dhpr 



1-S2 



dhpr 



(4) 
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^L =k sirsi dhpr -ksfrs2 dhpr (5) 



Expressions for the rate constants: 



CaM] dhp -100.0 



(6) 



fes^T =18.0 (7) 

Gating scheme for the 6-state Markovian model for the DHP-sensitive L-type Ca 2+ 
release channel: 



C^dhpr ~ 1 ~ (C^dhpr + 0 2 dhpr + °^dhpr) 

-iC4 dhpr ) 



(8) 



dCl dhpr _ h dhpr ni h dh P r C^ 

-(k?r+k^)c5 dhpr 



^ -M2 ^ L dhpr "•" "42 ^^dhpr 



-(fe 2 r+fe 2 T)02 dhpr 



dC4 dhpr _ h dhpr C r , udhpr m 

— - ^54 ^^dhpr + ^24 UZ dhpr ^) 

-(fe 4 r+fe 4 r)C4^ r 



dC6 ^=fe 3 ro3 dhpr -fe 6 rc6 dhpr d3) 



The open state 03 dhpr accounts for the increased tail current produced as the result 
of a large depolarization. Expressions for rate constants: 

fcff =1000.0 (14) 



feff = 600000.0 



(15) 
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h dhpr h dhpr h dhpr 
idhpr _ "45 K 52 ^24 
54 , dhfin dhfir 



k dhpr k dhpr 



(16) 



ydhpridhpridhpr 
fedhpr _ "51 K \2 K 25 

15 kf^ r kf^ r 



(17) 



Case 1: With CICR (control) 



k dhpr = 



300.3808- 



301.0817 

7 v-13.5918 ^ 
l + e { 8.85 J 



(18) 



k dhpr = 



3359.8754 + 



966.95 



l + ^l 1 - 45 



66 
.4585 



(19) 



§ = 550 + 6xCaMKII act 
+CaN nrf 



(20) 



336160 xS2 dhpr 



(21) 



5939.4 + 



306806.8 



1 + e 



v+1.8716 ^ 
1.3072 



(22) 



„„™^r- 0.48961 
0.02925 + - 



1 + e 



v+12.2249 
0.974 



(23) 



1350 1700 
1352 t - + 



1 + e 



1.8 



1 + 0 



(v+5) 



(24) 



1030.0575- 



713.1966 

1 + ^5.0) 



(25) 
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* 5 r= ,^ 



383.5435 0.16725 



8.3702 1 ( v-31.8252 



047 J 1 + £ V 0.01075 



26 

378.7085 

-373.015 



( v-7.8953 
l + 6 7691 



fe 3 t f = 



f \ 

loo- 1000 



( v-41.0 



(27) 



600.0 



v-41.0 
1.0 



(28) 



♦ Case 2: With Ryanodine applied (no RyR release) 



k dhpr = 



30392.87- 



30394.73 



83.5975 
28.5662 



f- 



(29) 



3568.74658 + 



21921.25344 



+ 24.9838 
7372 



( v+: 



(30) 



§ =550 + 6xCoA4KII flCt 
+CaN /irt 



(31) 



336160 XS2^ 



(32) 



0.0893 + 



26.3268 

v+15.0 ^ 
1.0 



1 + e 



(33) 



32.6642 



63.808 



1 + e 



-14.4897 
0.3131 



1 + gV 



v-TlAW 
0.0854 



(34) 
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d)ip _OT4152_ 
K 52 -Z.0J+ / v+27 945 

l + e { 1-46976 

2.43098 



f v-40.25 



(35) 



^=1908-™ 

15 I v+15.0 

l + el^ 99 * 
879.543 



1 + e 



v-11.1635 
0.3238 



(36) 



fe 3 T =145107.50934- 

f v-1037.5224 ^ 
1 + 312.91482 J 

5537.7875 



l + 02893 J 



(37) 



k dhpr = 



27.651 



11.488 



J v-59.48 "j ( v+15.1148 "j 

1 + e \ 4.444 J 1 + 0.9423 J 



1 + e 



t.444 

16. 3777 

v+15.01 
0.6351 



+ 5.3384 



(38) 



100 



100.0 



( v-41.0 

1 + el^^ 



(39) 



600.0 



f v-41.0 "j 
v l + ^l 10 

Case 3: With Ca 2+ substituted with Ba 2+ 
ktr= 3028.7604 6279 ; 833 



1 + e 



v-58.1484 
9.5348 



1 + e 



3259.706 



v-40.2856 
0.2272 



(40) 



(41) 
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k dhpr _ 26296.0826 + 4" ) 8.8 1 5 2 



( v+17.40398 _( v+1.4892 

1 + e { 0.22797 J 1 + 4.9581 

1544.9745 

-808.0319 



_r v-i2 
i + <A 3 

kfi r = (o.o) 



32.7305 



6234 



fc^= -5883.3476 + - 5973 3014 



( v-41.4221 ' 
1 + 05703 , 



528.171 



1 + 0 



v-20.4274 
0.187 



fe 5 d f =175.78742 ^Sr- 

b2 f v-32.3072 



1 + 0V 



^ v 4.4783 



229.8479 



v-36.9508 
5.4521 



l + e 

2601.9597 2647.5114 



123 f v-19.882 "I 

l + e { 0.4459 J 1 + g 



v+15.6872 
2.1947 



2501.6394 

+ 2647.2042 



l + e 



v-39.8404 
0.3926 



fe ^pr 1683.7686 2 1 6 1 >699 



_( v-20.0452 "| f v-360.2862 "| 

1 + e [ 0.40326 J 1 + ^ 1879.441 J 

5926.5225 „ _ 

+ 4129.8882 



1 + 0 



v+1.7842 
13.43183 



fe ijpr 



11.1993 36.2341 



_( v-35.2478 "| J v+23. 

1 + 0.2466 J 1 + 

1Z163 +37.7673 



9615 
8653 



1 + 0 



v-27.1529 
0.9455 
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k dhpr = 



f \ 

100.0 

100- 



( v-41.0 



(49) 



fc gpr = 



600.0 



1 + e 



v-41.0 
1.0 



(50) 



Uptake of Ca 2+ from the cytosol into the LSR 

Ca 2+ fluxes from cytosol to SERCA and SERCA to LSR 

J cyt,serca = ^ cyt, serca X V^ a \myo X SERCA m 
~k cyt, serca X [Cd ] m y 0 X [Cd ] serca 
~^ serca, cyt x ] serca 



(51) 



J serca, sr ~ ^sr, serca X + ]lSK X + \serca 
~^sr, serca X [C^^LSR X SERCA m 
+kcyt,sr X [C a \myo 



(52) 



r _ [ v ?FV 

1 cyt,serca J cyt,serca serca 



j — j x 9 ry 

1 serca, sr J serca, sr serca 

Differential equation for Ca 2+ buffered by the SERCA protein: 



d[Cd 



2+1 



dt 



' cyt, serca J serca, sr 



Expressions for the rate constants for Ca 2+ binding to/release from SERCA: 



Kcyt,serca£.07 + 5500 X CdMKII act ) 



v cyt, serca iV cyt, 



X - K ( prf wd ^ 

iV cyt, serca ~ LX cyt,serca\ ^^50 ) 



kserca.sr = ^serca.sM' 07 + 5500 X CaMKU act) 



(53) 
(54) 

(55) 

(56) 
(57) 
(58) 



_ serca, sr 
~ cn n bwd\2 
l^ C 50 J 



(59) 



Krishna et al. Theoretical Biology and Medical Modelling 2010, 7:43 
http://www.tbiomed.eom/content/7/1/43 



Page 53 of 66 



Affinities for the forward and backward Ca 2+ fluxes: 
ECl% d = 0.015(1 + PSRxPLB dp ) 



1.027 



1 + 5500 xCaMKII, 



(60) 



EC h 5 Q d = 1250 - 1 1 10 x PSR x PLB dp 



Differential equation for phospholamban 

dPI ^t = k™PLB p 



dt 



-V 2 f (5500CaMKII act + 
2.5CaN act +PKA act ) 2 PLB dp 



(61) 



(62) 



PLB p =l-PLB dp 



(63) 



Ca pump in SL 



IpMCA - IpMCA 



[Ca 2 



imyo 



kmpca + [Ca 2+ ] myo 



(64) 



Na + /Ca 2+ exchanger 



_ NUMl x (NUM2 - NUM3) 

* NaCa - ' — n 



( 



-0.65FV \ 



l + 0.27e 



(65) 



NUMl = R 



NaCa 



' K * 
A mAllo 



1 + 

V V 



[Ca 2+ ] NaCa 



J ) 



(66) 



0.35FV 

NUM2 = e «r {[Na + ] NaCa ) 3 [Ca 2+ ] 0 



(67) 



0.65FV 

NUMS = e «r {[Na + ] 0 ) 3 [Ca 2+ } NaCa 



(68) 
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P - K m caoi[Na + ] NaCa ) 3 + ^ Nao [Ca 2+ ] NaCa 

( r^.2+i \ 



+ K mNai[ Ca2+ ] 



1 + 



[Ca 2+ ] NaCa 



K 



mCai 



+K„q,,([N«*U 5 

+ (K]J'[ C » J *] 



( / A 3 



mNai 



(69) 



where [C<2 2+ ]7v aC «> [Na + ] NaCa are Concentrations at the mouth of the NaCa- 
exchanger. 

is obtained by dividing the above net chan- 



The corresponding unitary current i^aCa ls obtair 
nel current luaCa by the number of dyadic units 



Na + /Cs + pump 



^NaCs ~ R NaCs^NaCs 



100([Na + ] myo ) 



[Cs + ] 0 
[Cs + ] 0 + kmcs 

+i \ 15 



(70) 



(Wa + ] myo ) 15 +kmna 
0.65558 

J v+53.353 "| 
0.18445 + 2 V 15 - 58 J y 

The corresponding unitary current i^aCs is obtained by dividing the above net chan- 
nel current ImuCs by the number of dyadic units N^ ya ^ 

Background Na + current 



r ( 

v In 

F 



V 



[Na + ] 0 



[Na 4 



Jrayo 



(71) 



Ca 2+ transfer from LSR to a single jSR 



f [Ca 2+ \ L sR -[Ca 2+ ]j S R ^ 



2FV jSR 



(72) 



Ca 2+ release from a unit jSR into a single DCU 



iryr =J ryr x(2FxxA 2 A z ) 



(73) 
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where, 



° 2 ryri[ Ca2+ ]jSR -[ CaZ+ ]ryr) P ryr 



2+1 



' ryr 



ttA 2 A z 



(74) 



Gating scheme for the 4-state Markovian model for the RyR-sensitive SR Ca 2+ 
release channel: 



C4 — 1 (Cl^ r + 02 + C3™.) 



(75) 



dCl 



ryr _ 



= k%C4 Tyr +k$02 T 



dt -« -ryr -21— ryr 



L ryr 



(76) 



— = fe 12 Cl ryr + /e 32 03^ 



ryr 



(77) 



dC3 



ryr 



(78) 



fe£ =[Ca 2+ ] 



2+i2 

ryr 



0.05 -- 



3.7465xl0~ 2 
l + e ^3 



CaMKII M 
12000 



(79) 



.ryr 
"21 



698.56 



618.56 



[ca 2+ ] fyr -2052.9 

1 + 0 



1 + - 



6.0 



[Ca 



2+1 



\ryr 



1 



(80) 



fe{^ = 8.748 xlO~ 2 [Ca 2+ ] var 



(81) 



i ryr 
«41 



2.0 ^ 



CaMKU^ 
12000 



(82) 



_ .ryr. t, ryr _ i ryr. 



(83) 
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_ iryriryr _ yyr . 



(84) 



var = 



Al ls -0.7013 

10.0 xe 003 



(85) 



where [Ca 2+ ] ryr is the Ca 2+ concentration at the mouth of the RyR channel on the 
dyadic side. In the presence of caffeine (CF, concentration in fiM): 



02 



( \ 
-0.522 



ryr 



[CF]-410 

l + e 2649 



+ 0.503 



(86) 



Gating scheme for the 6-state Markovian model for the Luminal Calcium sensor. 

/4 is =l-(Al, s +i2 is +/3, 5 ) (87) 



B5 is =l-(/2, 5 +/3, s + B6, s ) 



^ = k^I4 u + k b 21 I2 b - (k[ s 2 B5 ls + k[\)Al k 



(88) 
(89) 



^ = k[ s 2 B5 ls Al ls + k k 2 B5 k I4 ls + k b 2 B ls 
-ik 2 \+k 2 s i+ k l l 3 + k 2 s 5 )I2 k 



(90) 



dt 



K,B5 h I4 h + k b 3 I2 k - (k k 4 + k b 2 ) I3 b 



(91) 



^ = kl 6 [Ca 2+ ] jSR B5 ls -k 1 i 5 B6 ls 

Ali s Fractional occupancy of RyR by Triadin/Junctin 

12 i s Fractional occupancy of RyR by Triadin/Junctin and Calsequestrin 

13 is Fractional occupancy of Triadin/Junctin by Calsequestrin 

74/ s Free Triadin/Junctin, B51s - Free Calsequestrin 

B6 is Fractional Occupancy of Calsequestrin by Calcium 

Ca 2+ dependent CaM mediated activation of CaMKII and CaN 

Ca 2+ binding to CaM and CaM buffering 



(92) 



CaM = CaM fl 



Ca 2 CaM -Ca 4 CaM - CaMB 



-Ca 2 CaMB -Ca 4 CaMB -CaMCaN 
-Ca 2 CaMCaN - Ca^CaMCaN 
-CaMKII tot ( Pl +P 3 + P 5 + P 4 ) 



(93) 
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B = B tot - [CaMB + Ca 2 CaMB + Ca 4 CaMB) 

R02 = ko2[Ca 2+ ] 2 x CaM - x Ca 2 CaM 

R24 = k 24 [Ca 2+ \ 2 x Ca 2 CaM - k 42 x Ca 4 CaM 

R02B = k% 2B [Ca 2+ ] 2 x CaMB - x Ca 2 CaMB 

R2AB = k™ B [Ca 2+ ] 2 x Ca 2 CaMB 
-k 42B x Ca 4 CaMB 

ROB = k^ on x CaA4 xB- k$f off x Ca 2 CaMB 
R2B = k 2Bon x Ca 2 CaM xB- k 2Boff x Ca 2 CaMB 
R4B = k 4Bon x Ca 4 CaM xB- k% Boff x Ca 4 CaMB 



dCa 2 CaM 



= R02 - R24 - R2B - R2CaN 



dt 



+CaMKLI m x [RCK 56 - RCK 21 ) 



dCa 4 CaM 
dt 



= R24 - R4CaN - RAB 



-CaMKII tot x [RCK 46 - RCK 23 ) 



dCa 4 CaM 
dt 



= R24 - R4CaN - R4B 



-CaMKII tot x {RCK 46 - RCK 23 ) 



dCaMB 
dt 



= ROB - R02B 



dCa 2 CaMB 
dt 



= R02B + R2B - R24B 



dCa 4 CaMB 
dt 



= R24B - RAB 



CaMKII activation: 
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T = P 3 +P 4 +P S + P 6 



(109) 



feff = 0.055 x T + 0.0074 xT J 



<34 



(no) 



+0.015XT' 



RCK M =k^xP 3 



K PP1 x PPl tot x P 4 
fe mppl +(CflMK// tot xP 4 ) 



(111) 



PCK 21 = ftff x Ca 2 CaM x P 2 



(112) 



KOC 13 = fcjf x [C<TT x Pj 

^32 ^ ^3 



(113) 



PCK 23 = fc 2 f x Ca 4 CaM x P 2 

^32 ^ "^3 



45 



M5 



-fc^f x[Ca^xP 5 



(114) 



(115) 



RCK 46 =k™xP 4 

-\ql x Ca 4 CaM x P 6 



PCK 56 = fc£f x P 5 

56 56 5 ^ 117 j 

-fc£f x Ca 2 CaM x P 6 



PCK 51 = feppiXm ^ XPs (118) 

fe mPP1 +(CaWII tot xP 5 ) 



PCK 62 = fepPiXPPl tQt xP 6 

fe mPP1+ (CaWII tot xP 6 ) 



dP 

= RCK 21 + PCK 51 - PCK 13 (120) 



dP 

— ^ = PCK 23 + PCK 13 - PCK 34 (121) 
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dP 

— ± = RCK 34 - RCK 46 - RCK 45 (122) 
dt 



dP 

= RCK 45 - RCK 56 - RCK 51 (123) 

dt 



dP 

— = RCK 46 - RCK 56 - RCK 62 (124) 
dt 



CaMKII act = 100 x (P 3 + P 4 + P 5 + P 6 ) (125) 
CaN activation: 



Ca 2 CaN = CaN tot - Ca 4 CaN - CaMCaN 
-Ca 2 CaMCaN - Ca 4 CaMCaN 



RCN Ca4 = kZn x [Ca 2+ ] 2 x Ca 2 CaN 
-kcaoff x Ca 4 CaN 



RCN 02 = k™ x [Ca 2+ ] 2 x CaMCaN 
-k 2 Q x Ca 2 CaMCaN 



RCN 24 = k 24 x [Ca 2+ ] 2 x Ca 2 CaMCaN 
-k 42 x Ca 4 CaMCaN 



RCN 0 = k^ n x CaM x Ca 4 CaN 



-feo^ x CaMCaN 



RCN 2 = k™ n x Ca 2 CaM x Ca 4 CaN 
-k™ ff x Ca 2 CaMCaN 



RCN 4 = k™ n x Ca 4 CaM x Ca 4 CaN 
-k™ ff x Ca 4 CaMCaN 



dt 



(126) 



(127) 



(128) 



(129) 



(130) 



(131) 



(132) 



dCa f aN = RCN Ca4 - RCN 0 - RCN 2 - RCN 4 (133) 



dCaMCaN = RCNq _ RCNq2 (134 ) 



(135) 
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dCa 4 CaMCaN ri ^ VT n o, T /iq£\ 
= RCN 4 + RCN 24 (136) 



CaN act = 100(Ca 4 CaMCaN + 0.1 Ca 2 CaMCaN) 

+ 10Ca 4 CaMCaN +0.lCaMCaN (137) 
+0.1Ca 4 CaN 



A2 - Differential equations for buffers used in the model 

Fluorescent indicator dye 

= kfluo3[Ca% 0 ]{[fluo3] tot - [CaF3]) 
-k- fluo3 [CaF3] 

Intracellular Ca 2+ buffering: 
calmodulin (bulkmyoplasm): 

^ = 200[Ca 2+ ] myo (l - OJ - 476.0 x O c (139) 

Troponin 



(138) 



(Fractional occupancy of troponin-Ca complex by Ca + ): 

= 78A[Ca 2+ } myo {l - O tc ) - 392.0 x O to (140) 
(Fractional occupancy of troponin-Mg complex by Ca 2+ ): 

fa = 200[Ca ] m y 0 (l — Otmgc ~ Otmgmg) (141) 

-6.6xO tmgc 

(Fractional occupancy of troponin-Mg complex by Mg 2 *): 

d °tmgmg = 2 Q[A4 2+i fl _ Q _ Q x 

^ ,u L iVI cS JrayoV^ w Lmgc w imgmg ) (142) 

-666.0 xO tmgmg 



A3 - Differential equations for ion concentrations used in the model 
Intracellular Ca 2+ concentration: 

1. Ca 2+ concentration in the cytosol 

d[Ca 2+ ] myo _ Ijyqd ~ IpMCA ~ Ica,L,SL 

dt ~ 2FV myo 

~Icyt,serca + ^ NaCa,SL (143) 

2FV 

myo 

dO dCaFl 
dt dt 
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where I dyad is the net integrated Ca 2+ flux diffusing out of all the dyadic units into 
the cytosol 



^ = 3.2^ + 6.4^^ + 1.8 
at at at 

2. Ca 2+ concentration in the jSR 
d[Ca 2+ 



dt 



dt 



jSR _ l tr l ryr _ 3 10 qq dB6 h 



2FV 



jSR 



dt 



3. Ca 2+ concentration in the LSR 



(144) 



(145) 



d[Cd ]isr _ ^cyt,serca 



' N dyadhr 



dt 2FV LSR 
4. Diffusion equation for Calcium in the dyadic space 



(146) 



d[Ca 



2+1 



\dyad 



dt 



:D Q ,V 2 [Ca 2+ ] dM + / f)T 



+ ^dhpr + J Naca + ^bnd 



(147) 



1 dhpr 



L Ca,L,TT 



N 



dyad 



5.1821X1CT 



(148) 



' NaCa.TT 



L NaCa,TT 



N 



dyad 



2x5.1821x10" 



(149) 



>bnd 



N h K h 



[K h +[Ca 2+ ] dyad ) 2 ^ 
{K l+ [Ca 2+ ] dyad ) 2 



d[Ca 2+ ] dyad 



dt 



(150) 



Intracellular Na + concentration: 

1. Na + concentration in the cytosol 



d[Na" 
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2. Na + concentration in the dyadic space 



d[Na + 
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dt 
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Intracellular Cs + concentration: 

1. Cs + concentration in the cytosol 
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2. Cs + concentration in the dyadic space 
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DHP-sensitive l Ca ,L channel; [Ca 2+ ] dyad : spatial Ca 2+ concentration in the dyad; [Ca 2+ ] jS R : luminal Ca 2+ concentration in 
the jSR; [Ca 2+ ] LS R : Ca 2+ concentration in the LSR; [Ca 2+ ] myo : cytosolic Ca 2+ concentration; [Ca 2+ ] serca : Ca 2+ 
concentration buffered by the SERCA protein; [Ca 2+ ] Q : extracellular Ca 2+ concentration; [Ca 2+ ] ryr : Ca 2+ concentration 
at the "mouth" of the RyR channel on the dyadic side; CaCM: Ca 2+ bound calmodulin; CaF3: Ca 2+ bound to the buffer 
Fluo3; CaM: calmodulin; CaMKII: Ca 2+ /calmodulin-dependent protein kinase II; CaMKII act : activated Ca 2+ /calmodu In- 
dependent protein kinase II; CaN: calcineurin; C a N act : activated calcineurin; CDF: calcium dependent facilitation; CDI: 
calcium dependent inactivation; CF: caffeine; CICR: calcium-induced calcium-release; CS: calsequestrin; [Cs + ] myo : Cs + 
concentration in the cytosol; [Cs + ] dyad : Cs + concentration in the dyadic space; DCU: dyadic coupling unit; DHP: 
dihydropyridine; DHPR: dihydropyridine receptor; E-C: excitation contraction; EC50bwd: affinity of backward Ca 2+ flux 
from LSR to cytosol; EC50fwd: affinity of forward Ca 2+ flux from cytosol to LSR; l Ca ,L ■ L-type Ca 2+ current; l Ca ,L,sL ■ 
sarcolemmal component of the l Ca ,L channel current; lca,L,n'- kai channel facing the dyadic space; l Cs : cesium current 



through the l Ca ,L channel; / c> 



, : Co uptake current directed from the cytosol to the SERCA; l Na : sodium current 



through the l Ca ,L channel; l Nab : background sodium current; l NaCa '■ sodium calcium exchanger current; l Na ca,sL '■ 
sarcolemmal component of the l NaCa exchanger current; lNaCa,rr- InoCo exchanger facing the dyadic space; l NaCs : 
sodium cesium pump current; l PMCA : plasma membrane Ca 2+ ATPase pump current; l ryr : Ca 2+ current due to CICR 



from an individual jSR; / 



: Co uptake current directed from the SERCA to the LSR; l tr : Ca current due to 



concentration gradient driven Ca transport from LSR to jSR; jSR: junctional portion of the sarcoplasmic reticulum; k d : 
dissociation constant; k mp : half saturation constant for the sarcolemmal Ca 2+ pump; LCC: L-type DHP-sensitive Ca 2+ 
channel; L-type: long lasting type; LSR: longitudinal portion of the sarcoplasmic reticulum; mM: milli molar; mV: milli 
volt; [Na + ] dyad : Na + concentration in the dyadic space; [Na + ] myo : Na + concentration in the cytosol; nM: nano molar; 
O c : fractional occupancy of calmodulin by Ca 2+ in the cytosol; O co/se : fractional occupancy of calsequestrin by Ca 2+ in 
the jSR; O tc : fractional occupancy of troponin-Ca sites by Ca 2+ in the cytosol; O tmgc : fractional ocupancy of troponin- 
Mg sites by Ca 2+ in the cytosol; O tmgmg : fractional occupancy of troponin-Mg sites by Mg 2+ in the cytosol; pC: pico 
coulomb; P Q : Open probability; PKA: protein kinase A; PLB dp : Unphosphorylated phospholamban; PLB P : 
phosphorylated phospholamban; PSR: phospholamban to SERCA ratio; Ry: ryanodine; RyR: ryanodine receptor; SERCA: 
sarcoplasmic reticulum Ca 2+ ATPase; SL: sarcolemma; SR: sarcoplasmic reticulum; Tc: Ca 2+ -specific troponin binding 
site; T-tubule: transverse tubules; TT: transverse tubules; VC: voltage clamp; VDI: voltage dependent inactivation; 
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